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1. 


INTRODUCTION 


The  Summer  Research  Program  (SRP),  sponsored  by  the  Air  Force  Office  of  Scientific  Research 
(AFOSR),  offers  paid  opportunities  for  university  faculty,  graduate  students,  and  high  school 
students  to  conduct  research  in  U.S.  Air  Force  research  laboratories  nationwide  during  the  summer. 

Introduced  by  AFOSR  in  1978,  this  innovative  program  is  based  on  the  concept  of  teaming 
academic  researchers  with  Air  Force  scientists  in  the  same  disciplines  using  laboratory  facilities  and 
equipment  not  often  available  at  associates'  institutions. 

The  Summer  Faculty  Research  Program  (SFRP)  is  open  annually  to  approximately  150  faculty 
members  with  at  least  two  years  of  teaching  and/or  research  experience  in  accredited  U.S.  colleges, 
universities,  or  technical  institutions.  SFRP  associates  must  be  either  U.S.  citizens  or  permanent 
residents. 

The  Graduate  Student  Research  Program  (GSRP)  is  open  annually  to  approximately  100  graduate 
students  holding  a  bachelor's  or  a  master's  degree;  GSRP  associates  must  be  U.S.  citizens  enrolled 
full  time  at  an  accredited  institution. 

The  High  School  Apprentice  Program  (HSAP)  annually  selects  about  125  high  school  students 
located  within  a  twenty  mile  commuting  distance  of  participating  Air  Force  laboratories. 

AFOSR  also  offers  its  research  associates  an  opportunity,  under  the  Summer  Research  Extension 
Program  (SREP),  to  continue  their  AFOSR-sponsored  research  at  their  home  institutions  through  the 
award  of  research  grants.  In  1994  the  maximum  amount  of  each  grant  was  increased  from  $20,000 
to  $25,000,  and  the  number  of  AFOSR-sponsored  grants  decreased  from  75  to  60.  A  separate 
annual  report  is  compiled  on  the  SREP. 

The  numbers  of  projected  summer  research  participants  in  each  of  the  three  categories  and  SREP 
“grants”  are  usually  increased  through  direct  sponsorship  by  participating  laboratories. 

AFOSR' s  SRP  has  well  served  its  objectives  of  building  critical  links  between  Air  Force  research 
laboratories  and  the  academic  community,  opening  avenues  of  communications  and  forging  new 
research  relationships  between  Air  Force  and  academic  technical  experts  in  areas  of  national 
interest,  and  strengthening  the  nation's  efforts  to  sustain  careers  in  science  and  engineering.  The 
success  of  the  SRP  can  be  gauged  from  its  growth  from  inception  (see  Table  1)  and  from  the 
favorable  responses  the  1997  participants  expressed  in  end-of-tour  SRP  evaluations  (Appendix  B). 

AFOSR  contracts  for  administration  of  the  SRP  by  civilian  contractors.  The  contract  was  first 
awarded  to  Research  &  Development  Laboratories  (RDL)  in  September  1990.  After  completion  of 
the  1990  contract,  RDL  (in  1993)  won  the  recompetition  for  the  basic  year  and  four  1-year  options. 


1 


2. 


PARTICIPATION  IN  THE  SUMMER  RESEARCH  PROGRAM 


The  SRP  began  with  faculty  associates  in  1979;  graduate  students  were  added  in  1982  and  high 
school  students  in  1986.  The  following  table  shows  the  number  of  associates  in  the  program  each 
year. 


1998 


85 


40 


88 


213 


Beginning  in  1993,  due  to  budget  cuts,  some  of  the  laboratories  weren’t  able  to  afford  to  fund  as 
many  associates  as  in  previous  years.  Since  then,  the  number  of  funded  positions  has  remained 
fairly  constant  at  a  slightly  lower  level. 


3.  RECRUITING  AND  SELECTION 

The  SRP  is  conducted  on  a  nationally  advertised  and  competitive-selection  basis.  The  advertising 
for  faculty  and  graduate  students  consisted  primarily  of  the  mailing  of  8,000  52-page  SRP  brochures 
to  chairpersons  of  departments  relevant  to  AFOSR  research  and  to  administrators  of  grants  in 
accredited  universities,  colleges,  and  technical  institutions.  Historically  Black  Colleges  and 
Universities  (HBCUs)  and  Minority  Institutions  (Mis)  were  included.  Brochures  also  went  to  all 
participating  USAF  laboratories,  the  previous  year's  participants,  and  numerous  individual 
requesters  (over  1000  annually). 

RDL  placed  advertisements  in  the  following  publications:  Black  Issues  in  Higher  Education,  Winds 
of  Change,  and  IEEE  Spectrum.  Because  no  participants  list  either  Physics  Today  or  Chemical  & 
Engineering  News  as  being  their  source  of  learning  about  the  program  for  the  past  several  years, 
advertisements  in  these  magazines  were  dropped,  and  the  funds  were  used  to  cover  increases  in 
brochure  printing  costs. 

High  school  applicants  can  participate  only  in  laboratories  located  no  more  than  20  miles  from  their 
residence.  Tailored  brochures  on  the  HSAP  were  sent  to  the  head  counselors  of  180  high  schools  in 
the  vicinity  of  participating  laboratories,  with  instructions  for  publicizing  the  program  in  their 
schools.  High  school  students  selected  to  serve  at  Wright  Laboratory's  Armament  Directorate 
(Eglin  Air  Force  Base,  Florida)  serve  eleven  weeks  as  opposed  to  the  eight  weeks  normally  worked 
by  high  school  students  at  all  other  participating  laboratories. 

Each  SFRP  or  GSRP  applicant  is  given  a  first,  second,  and  third  choice  of  laboratory.  High  school 
students  who  have  more  than  one  laboratory  or  directorate  near  their  homes  are  also  given  first, 
second,  and  third  choices. 

Laboratories  make  their  selections  and  prioritize  their  nominees.  AFOSR  then  determines  the 
number  to  be  funded  at  each  laboratory  and  approves  laboratories'  selections. 

Subsequently,  laboratories  use  their  own  funds  to  sponsor  additional  candidates.  Some  selectees  do 
not  accept  the  appointment,  so  alternate  candidates  are  chosen.  This  multi-step  selection  procedure 
results  in  some  candidates  being  notified  of  their  acceptance  after  scheduled  deadlines.  The  total 
applicants  and  participants  for  1998  are  shown  in  this  table. 
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i  1998  Applicants  and  Participants 

PARTICIPANT 

TOTAL 

SELECTEES 

DECLINING 

CATEGORY 

APPLICANTS 

SELECTEES 

SFRP 

382 

85 

13 

(HBCU/MI) 

(0) 

(0) 

(Q) _ 

GSRP 

130 

40 

7 

(HBCU/MI) 

(0) 

(Q) 

(0) 

HSAP 

328 

88 

22 

TOTAL 

840 

213 

42 

4.  SITE  VISITS 

During  June  and  July  of  1998,  representatives  of  both  AFOSR/NI  and  RDL  visited  each 
participating  laboratory  to  provide  briefings,  answer  questions,  and  resolve  problems  for  both 
laboratory  personnel  and  participants.  The  objective  was  to  ensure  that  the  SRP  would  be  as 
constructive  as  possible  for  all  participants.  Both  SRP  participants  and  RDL  representatives  found 
these  visits  beneficial.  At  many  of  the  laboratories,  this  was  the  only  opportunity  for  all  participants 
to  meet  at  one  time  to  share  their  experiences  and  exchange  ideas. 


5.  HISTORICALLY  BLACK  COLLEGES  AND  UNIVERSITIES  AND  MINORITY 
INSTITUTIONS  (HBCU/MIs) 

Before  1993,  an  RDL  program  representative  visited  from  seven  to  ten  different  HBCU/MIs 
annually  to  promote  interest  in  the  SRP  among  the  faculty  and  graduate  students.  These  efforts  were 
marginally  effective,  yielding  a  doubling  of  HBCI/MI  applicants.  In  an  effort  to  achieve  AFOSR  s 
goal  of  10%  of  all  applicants  and  selectees  being  HBCU/MI  qualified,  the  RDL  team  decided  to  try 
other  avenues  of  approach  to  increase  the  number  of  qualified  applicants.  Through  the  combined 
efforts  of  the  AFOSR  Program  Office  at  Bolling  AFB  and  RDL,  two  very  active  minority  groups 
were  found,  HACU  (Hispanic  American  Colleges  and  Universities)  and  AISES  (American  Indian 
Science  and  Engineering  Society).  RDL  is  in  communication  with  representatives  of  each  of  these 
organizations  on  a  monthly  basis  to  keep  up  with  the  their  activities  and  special  events.  Both 
organizations  have  widely-distributed  magazines/quarterlies  in  which  RDL  placed  ads. 

Since  1994  the  number  of  both  SFRP  and  GSRP  HBCU/MI  applicants  and  participants  has  increased 
ten-fold,  from  about  two  dozen  SFRP  applicants  and  a  half  dozen  selectees  to  over  100  applicants 
and  two  dozen  selectees,  and  a  halt-dozen  GSRP  applicants  and  two  or  three  selectees  to  18 
applicants  and  7  or  8  selectees.  Since  1993,  the  SFRP  had  a  two-fold  applicant  increase  and  a  two¬ 
fold  selectee  increase.  Since  1993.  the  GSRP  had  a  three-fold  applicant  increase  and  a  three  to  four¬ 
fold  increase  in  selectees. 
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In  addition  to  RDL's  special  recruiting  efforts,  AFOSR  attempts  each  year  to  obtain  additional 
funding  or  use  leftover  funding  from  cancellations  the  past  year  to  fund  HBCU/MI  associates. 


SRP  HBCU/MI  Participation,  By  Year 

YEAR 

SFRP 

GSRP 

Applicants 

Participants 

Applicants 

Participants 

1985 

76 

23 

15 

11 

1986 

70 

18 

20 

10 

1987 

82 

32 

32 

10 

1988 

53 

17 

23 

14 

1989 

39 

15 

13 

4 

1990 

43 

14 

17 

3 

1991 

42 

13 

8 

5 

1992 

70 

13 

9 

5 

1993 

60 

13 

6 

2 

1994 

90 

16 

11 

6 

1995 

90 

21 

20 

8 

1996 

119 

27 

18 

7 

6.  SRP  FUNDING  SOURCES 

Funding  sources  for  the  1998  SRP  were  the  AFOSR-provided  slots  for  the  basic  contract  and 
laboratory  funds.  Funding  sources  by  category  for  the  1998  SRP  selected  participants  are  shown 
here. 
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1998  SRP  FUNDING  CATEGORY 

SFRP 

GSRP 

HSAP 

AFOSR  Basic  Allocation  Funds 

67 

38 

75 

USAF  Laboratory  Funds 

17 

2 

13 

Slots  Added  by  AFOSR 

0 

0 

0 

(Leftover  Funds) 

HBCU/MI  By  AFOSR 

0 

0 

N/A 

(Using  Procured  Addn’l  Funds) 

TOTAL 

84 

40 

88 

7.  COMPENSATION  FOR  PARTICIPANTS 

Compensation  for  SRP  participants,  per  five-day  work  week,  is  shown  in  this  table. 


1998  SRP  Associate  Compensation 


PARTICIPANT  CATEGORY 

1991 

1992 

1993 

1994 

1995 

1996 

1997 

1998 

Faculty  Members 

$690 

$718 

$740 

$740 

$770 

$793 

Graduate  Student 
(Master's  Degree) 

$425 

$442 

$455 

$455 

$455 

$470 

$470 

$484 

Graduate  Student 

(Bachelor's  Degree) 

$365 

$380 

$391 

$391 

$391 

$400 

$400 

$412 

High  School  Student 
(First  Year) 

$200 

$200 

$200 

$200 

$200 

High  School  Student 
(Subsequent  Years) 

$240 

$240 

$240 

$240 

$240 

/ 

The  program  also  offered  associates  whose  homes  were  more  than  50  miles  from  the  laboratory  an 
expense  allowance  (seven  days  per  week)  of  $52/day  for  faculty  and  $41/day  for  graduate  students. 
Transportation  to  the  laboratory  at  the  beginning  of  their  tour  and  back  to  their  home  destinations  at 
the  end  was  also  reimbursed  for  these  participants.  Of  the  combined  SFRP  and  GSRP  associates, 

65  %  claimed  travel  reimbursements  at  an  average  round-trip  cost  of  $730. 

Faculty  members  were  encouraged  to  visit  their  laboratories  before  their  summer  tour  began.  All 
costs  of  these  orientation  visits  were  reimbursed.  Forty-three  percent  (85  out  of  188)  of  faculty 
associates  took  orientation  trips  at  an  average  cost  of  $449.  By  contrast,  in  1993,  58  %  of  SFRP 
associates  elected  to  take  an  orientation  visits  at  an  average  cost  of  $685;  that  was  the  highest 
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percentage  of  associates  opting  to  take  an  orientation  trip  since  RDL  has  administered  the  SRP,  and 
the  highest  average  cost  of  an  orientation  trip. 

Program  participants  submitted  biweekly  vouchers  countersigned  by  their  laboratory  research  focal 
point,  and  RDL  issued  paychecks  so  as  to  arrive  in  associates'  hands  two  weeks  later. 

This  is  the  third  year  of  using  direct  deposit  for  the  SFRP  and  GSRP  associates.  The  process  went 
much  more  smoothly  with  respect  to  obtaining  required  information  from  the  associates,  about  15% 
of  the  associates’  information  needed  clarification  in  order  for  direct  deposit  to  properly  function  as 
opposed  to  7%  from  last  year.  The  remaining  associates  received  their  stipend  and  expense 
payments  via  checks  sent  in  the  US  mail. 

HSAP  program  participants  were  considered  actual  RDL  employees,  and  their  respective  state  and 
federal  income  tax  and  Social  Security  were  withheld  from  their  paychecks.  By  the  nature  of  their 
independent  research,  SFRP  and  GSRP  program  participants  were  considered  to  be  consultants  or 
independent  contractors.  As  such,  SFRP  and  GSRP  associates  were  responsible  for  their  own 
1  income  taxes,  Social  Security,  and  insurance. 

8.  CONTENTS  OF  THE  1998  REPORT 

The  complete  set  of  reports  for  the  1998  SRP  includes  this  program  management  report  (Volume  1) 
augmented  by  fifteen  volumes  of  final  research  reports  by  the  1998  associates,  as  indicated  below: 


1998  SRP  Final  Report  Volume  Assignments 


LABORATORY 

SFRP 

GSRP 

HSAP 

Armstrong 

2 

7 

12 

Phillips 

3 

8 

13 

Rome 

4 

9 

14 

Wright 

5A,  5B 

10 

15 

AEDC,  ALCs,  USAFA,  WHMC 

6 

11 

7 


APPENDIX  A  ~  PROGRAM  STATISTICAL  SUMMARY 


A.  Colleges/Universities  Represented 

Selected  SFRP  associates  represented  169  different  colleges,  universities,  and  institutions, 
GSRP  associates  represented  95  different  colleges,  universities,  and  institutions. 

B.  States  Represented 

SFRP  -Applicants  came  from  47  states  plus  Washington  D.C.  Selectees  represent  44  states. 
GSR?  -  Applicants  came  from  44  states.  Selectees  represent  32  states. 

HSAP  -  Applicants  came  from  thirteen  states.  Selectees  represent  nine  states. 


Total  Number  of  Participants 

SFRP 

85 

GSRP 

40 

HSAP 

88 

TOTAL 

213 

Degrees  Represented 

SFRP 

GSRP 

TOTAL 

Doctoral 

83 

0 

83 

Master's 

1 

3 

4 

Bachelor's 

0 

22 

22 

TOTAL 

186 

25 

109 

A-l 


SFRP  Academic  Titles 


Assistant  Professor 

36 

Associate  Professor 

34 

Professor 

15 

Instructor 

0 

Chairman 

0 

Visiting  Professor 

0 

Visiting  Assoc.  Prof. 

0 

Research  Associate 

0 

TOTAL 

85 

Source  of  Learning  About  the  SRP 


Category 

Applicants 

Selectees 

Applied/participated  in  prior  years 

177 

47 

Colleague  familiar  with  SRP 

104 

24 

Brochure  mailed  to  institution 

101 

21 

Contact  with  Air  Force  laboratory 

101 

39 

IEEE  Spectrum 

12 

1 

BIIHE 

4 

0 

Other  source 

117 

30 

TOTAL 


616 


162 


APPENDIX  B  -  SRP  EVALUATION  RESPONSES 


1.  OVERVIEW 

Evaluations  were  completed  and  returned  to  RDL  by  four  groups  at  the  completion  of  the  SRP.  The 
number  of  respondents  in  each  group  is  shown  below. 


Table  B-l.  Total  SRP  Evaluations  Received 


Evaluation  Group 

Responses 

SFRP  &  GSRPs 

100 

HSAPs 

75 

USAF  Laboratory  Focal  Points 

84 

USAF  Laboratory  HSAP  Mentors 

6 

All  groups  indicate  unanimous  enthusiasm  for  the  SRP  experience. 


The  summarized  recommendations  for  program  improvement  from  both  associates  and  laboratory 
personnel  are  listed  below: 


A.  Better  preparation  on  the  labs’  part  prior  to  associates'  arrival  (i.e.,  office  space, 
computer  assets,  clearly  defined  scope  of  work). 

B.  Faculty  Associates  suggest  higher  stipends  for  SFRP  associates. 

C.  Both  HSAP  Air  Force  laboratory  mentors  and  associates  would  like  the  summer  tour 

extended  from  the  current  8  weeks  to  either  10  or  11  weeks;  the  groups  state  it  takes 
4-6  weeks  just  to  get  high  school  students  up-to-speed  on  what’s  going  on  at 

laboratory.  (Note:  this  same  argument  was  used  to  raise  the  faculty  and  graduate 

student  participation  time  a  few  years  ago.) 
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2.  1998  USAF  LABORATORY  FOCAL  POINT  (LFP)  EVALUATION  RESPONSES 


The  summarized  results  listed  below  are  from  the  84  LFP  evaluations  received. 
1.  LFP  evaluations  received  and  associate  preferences: 


Table  B-2.  Air  Force  LFP  Evaluation  Responses  (By  Type) 


How  Many  Associates  Would  You  Prefer  To  Get 

_ r  ' 

7 

(%  Response) 

SFRP 

GSRP  (w/Univ  Professor) 

GSRP  (w/o  Univ  Professor) 

Lab 

Evals 

Reev’d 

0 

1 

2 

3+ 

0 

1 

2 

3+ 

0 

i 

2 

3+ 

AEDC 

0 

- 

- 

- 

- 

- 

- 

- 

- 

WHMC 

0 

- 

- 

- 

- 

- 

" 

- 

- 

AL 

7 

28 

28 

28 

14 

54 

14 

28 

0 

86 

0 

14 

0 

USAFA 

1 

0 

100 

0 

0 

100 

0 

0 

0 

0 

100 

0 

0 

PL 

25 

40 

40 

16 

4 

88 

12 

0 

0 

84 

12 

4 

0 

RL 

5 

60 

40 

0 

0 

80 

10 

0 

0 

100 

0 

0 

0 

WL 

46 

30 

43 

20 

6 

78 

17 

4 

0 

93 

4 

2 

0 

Total 

84 

32% 

50% 

13% 

5% 

80% 

11% 

6% 

0% 

73% 

23% 

4% 

0% 

LFP  Evaluation  Summary.  The  summarized  responses,  by  laboratory,  are  listed  on  the  following 
page.  LFPs  were  asked  to  rate  the  following  questions  on  a  scale  from  1  (below  average)  to  5 
(above  average). 

2.  LFPs  involved  in  SRP  associate  application  evaluation  process: 

a.  Time  available  for  evaluation  of  applications: 

b.  Adequacy  of  applications  for  selection  process: 

3.  Value  of  orientation  trips: 

4.  Length  of  research  tour: 

5  a.  Benefits  of  associate's  work  to  laboratory: 
b.  Benefits  of  associate's  work  to  Air  Force: 

6.  a.  Enhancement  of  research  qualifications  for  LFP  and  staff: 

b.  Enhancement  of  research  qualifications  for  SFRP  associate: 

c.  Enhancement  of  research  qualifications  for  GSRP  associate: 

7.  a.  Enhancement  of  knowledge  for  LFP  and  staff: 

b.  Enhancement  of  knowledge  for  SFRP  associate: 

c.  Enhancement  of  knowledge  for  GSRP  associate: 

8.  Value  of  Air  Force  and  university  links: 

9.  Potential  for  future  collaboration: 

10.  a.  Your  working  relationship  with  SFRP: 
b.  Your  working  relationship  with  GSRP: 

11.  Expenditure  of  your  time  worthwhile: 

(Continued  on  next  page) 
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12.  Quality  of  program  literature  for  associate: 

13.  a.  Quality  of  RDL’s  communications  with  you: 

b.  Quality  of  RDL's  communications  with  associates: 

14.  Overall  assessment  of  SRP: 


Table  B 

!-3.  Laboratory  Focal  Point  Reponses  to  above  questions 

AEDC 

AL 

USAF 

A 

PL 

RL 

WHMC 

WL 

it  Evals  Recv  ’d 

7 

1 

14 

5 

0 

46 

Question  if 

2 

- 

86  % 

0  % 

88  % 

80  % 

- 

85  % 

2a 

- 

4.3 

n/a 

3.8 

4.0 

- 

3.6 

2b 

- 

4.0 

n/a 

3.9 

4.5 

- 

4.1 

3 

- 

4.5 

n/a 

4.3 

4.3 

- 

3.7 

4 

- 

4.1 

4.0 

4.1 

4.2 

- 

3.9 

5a 

- 

4.3 

5.0 

4.3 

4.6 

- 

4.4 

5b 

- 

4.5 

n/a 

4.2 

4.6 

- 

4.3 

6a 

- 

4.5 

5.0 

4.0 

4.4 

- 

4.3 

6b 

_ 

4.3 

n/a 

4.1 

5.0 

- 

4.4 

6c 

- 

3.7 

5.0 

3.5 

5.0 

- 

4.3 

7a 

- 

4.7 

5.0 

4.0 

4.4 

- 

4.3 

7b 

- 

4.3 

n/a 

4.2 

5.0 

- 

4.4 

7c 

- 

4.0 

5.0 

3.9 

5.0 

- 

4.3 

8 

- 

4.6 

4.0 

4.5 

4.6 

- 

4.3 

9 

- 

4.9 

5.0 

4.4 

4.8 

- 

4.2 

10a 

- 

5.0 

n/a 

4.6 

4.6 

- 

4.6 

10b 

- 

4.7 

5.0 

3.9 

5.0 

- 

4.4 

11 

- 

4.6 

5.0 

4.4 

4.8 

- 

4.4 

12 

- 

4.0 

4.0 

4.0 

4.2 

- 

3.8 

13a 

- 

3.2 

4.0 

3.5 

3.8 

- 

3.4 

13b 

- 

3.4 

4.0 

3.6 

4.5 

- 

3.6 

14 

- 

4.4 

5.0 

4.4 

4.8 

- 

4.4 
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3.  1998  SFRP  &  GSRP  EVALUATION  RESPONSES 


The  summarized  results  listed  below  are  from  the  120  SFRP/GSRP  evaluations  received. 

Associates  were  asked  to  rate  the  following  questions  on  a  scale  from  1  (below  average)  to  5  (above 
average)  -  by  Air  Force  base  results  and  over-all  results  of  the  1998  evaluations  are  listed  after  the 
questions. 

1.  The  match  between  the  laboratories  research  and  your  field: 

2.  Your  working  relationship  with  your  LFP: 

3.  Enhancement  of  your  academic  qualifications: 

4.  Enhancement  of  your  research  qualifications: 

5.  Lab  readiness  for  you:  LFP,  task,  plan: 

6.  Lab  readiness  for  you:  equipment,  supplies,  facilities: 

7.  Lab  resources: 

8.  Lab  research  and  administrative  support: 

9.  Adequacy  of  brochure  and  associate  handbook: 

10.  RDL  communications  with  you: 

11.  Overall  payment  procedures: 

12.  Overall  assessment  of  the  SRP: 

13.  a.  Would  you  apply  again? 

b.  Will  you  continue  this  or  related  research? 

14.  Was  length  of  your  tour  satisfactory? 

15.  Percentage  of  associates  who  experienced  difficulties  in  finding  housing: 

16.  Where  did  you  stay  during  your  SRP  tour? 

a.  At  Home: 

b.  With  Friend: 

c.  On  Local  Economy: 

d.  Base  Quarters: 

17.  Value  of  orientation  visit: 

a.  Essential: 

b.  Convenient: 

c.  Not  Worth  Cost: 

d.  Not  Used: 

SFRP  and  GSRP  associate’s  responses  are  listed  in  tabular  format  on  the  following  page. 
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Table  B-4.  1997  SFRP  &  GSRP  Associate  Responses  to  SRP  Evaluation 


4.  1998  USAF  LABORATORY  HSAP  MENTOR  EVALUATION  RESPONSES 

Not  enough  evaluations  received  (5  total)  from  Mentors  to  do  useful  summary . 
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5.  1998  HSAP  EVALUATION  RESPONSES 

The  summarized  results  listed  below  are  from  the  23  HSAP  evaluations  received. 


HSAP  apprentices  were  asked  to  rate  the  following  questions  on  a  scale  from 
1  (below  average)  to  5  (above  average) 

1.  Your  influence  on  selection  of  topic/type  of  work. 

2.  Working  relationship  with  mentor,  other  lab  scientists. 

3.  Enhancement  of  your  academic  qualifications. 

4.  Technically  challenging  work. 

5.  Lab  readiness  for  you:  mentor,  task,  work  plan,  equipment. 

6.  Influence  on  your  career. 

7.  Increased  interest  in  math/science. 

8.  Lab  research  &  administrative  support. 

9.  Adequacy  of  RDL’s  Apprentice  Handbook  and  administrative  materials. 

10.  Responsiveness  of  RDL  communications. 

11.  Overall  payment  procedures. 

12.  Overall  assessment  of  SRP  value  to  you. 

13.  Would  you  apply  again  next  year?  Yes  (92  %) 

14.  Will  you  pursue  future  studies  related  to  this  research?  Yes  (68  %) 

15.  Was  Tour  length  satisfactory?  Yes  (82%) 


□ 

Arnold 

Brooks 

msg 

■33m 

Griffiss 

Hanscom 

Tyndall 

WPAFB 

Totals 

mm 

5 

19 

7 

15 

13 

2 

7 

5 

40 

113 

i 

2.8 

3.3 

3.4 

3.5 

3.4 

3.2 

3.6 

3.6 

3.4 

2 

4.4 

4.6 

4.5 

4.8 

4.6 

4.4 

4.0 

4.6 

4.6 

4.0 

4.2 

4.1 

4.3 

4.5 

mm 

4.3 

4.6 

4.4 

4.4 

3.6 

3.9 

4.0 

4.5 

4.2 

Bl 

4.6 

3.8 

4.3 

4.2 

5 

4.4 

4.1 

4.1 

3.9 

3.6 

4.0 

6 

3.2 

3.6 

|pii 

H9 

3.8 

mSM 

3.3 

3.8 

■M 

7 

2.8 

4.1 

3.9 

3.9 

5.0 

3.6 

4.0 

■El 

3.9 

8 

3.8 

4.1 

4.3 

4.0 

4.0 

4.3 

3.8 

m 

4.2 

9 

4.4 

3.6 

4.1 

4.1 

3.9 

■El 

3.7 

3.8 

10 

3.8 

4.1 

3.7 

3.9 

eh 

3.8 

3.8 

11 

4.2 

4.2 

3.8 

■9 

3.7 

2.6 

12 

4.5 

4.6 

Hi 

4.6 

4.2 

Numbers  below  are  percentages 

13 

60% 

95% 

100% 

100 

85% 

100% 

100% 

100% 

90% 

92% 

% 

14 

20% 

80% 

71% 

80% 
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A  STUDY  OF  OPTIMAL  FINITE-THRUST  SPACECRAFT 
TRAJECTORIES  FOR  THE  TECHS  AT  21  MISSION 


Benjamin  J.  Bernocco 
Graduate  Student 

Department  of  Aerospace  Engineering 
The  Pennsylvania  State  University 

Abstract 

A  preliminary  study  of  optimal  finite-thrust  spacecraft  trajectories  for  the  TechSat  21  mission  was 
examined.  The  study  looked  only  at  one  satellite  orbiting  in  the  ecliptic  plane,  i.e.  two-dimensional  motion,  with  no 
explicit  orbit  perturbations  included.  Optimization  of  the  spacecraft  trajectory  was  accomplished  using  a  direct 
collocation  with  nonlinear  programming  approach.  This  approach  requires  that  the  equations  of  motion  and 
integrals  of  the  motion  for  the  spacecraft  trajectory  be  calculated  from  state  vector  values  and  then  feed  into  an 
optimization  subroutine.  The  optimization  subroutine  changes  the  state  vector  values  to  minimize  the  objective 
function,  i.e.  the  thruster  firing  time,  which  is  subject  to  constraints  dictated  by  the  mission  requirements.  An 
example  simulation  for  this  mission  was  performed.  The  trajectory  was  divided  up  into  a  sequence  of  twelve  thrust 
arcs  separated  by  twelve  coast  arcs.  Constraints  were  placed  on  the  objective  function  to  insure  that  the  state  vector 
values  yielded  results  that  represented  an  actual,  predetermined  trajectory.  A  general  perturbation  was  incorporated 
into  the  problem  by  including  thrust  acceleration  terms  in  the  equations  of  motion.  The  initial  conditions  used  to 
start  the  process  were  chosen  to  be  the  state  vector  values  for  a  pure  Keplerian  orbit  with  no  perturbations.  After 
only  two  iterations  the  code  converged  to  an  optimized  solution.  This  solution  was  very  close  to  the  initial 
condition  with  only  the  thrust  level  and  thruster  firing  times  being  altered  by  the  program.  It  is  a  reasonable  solution 
since  the  thrusters  are  only  used  to  maintain  the  same  orbit,  hence  most  of  the  state  vector  values  would  not  change. 
The  solution  was  plugged  back  into  the  code  as  initial  conditions  and  the  code  yielded  the  same  solution  after  one 
iteration.  Although  this  does  not  guarantee  a  global  minimum  was  found,  it  instills  confidence  that  the  solution  is  a 
valid  one.  Future  work  to  be  done  includes  expanding  the  code  to  three  dimensions,  including  perturbation  effects 
like  atmospheric  drag,  luni-solar,  and  the  Earth’s  oblateness,  and  including  multiple  satellites  in  the  analysis. 
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A  STUDY  OF  OPTIMAL  FINITE-THRUST  SPACECRAFT 
TRAJECTORIES  FOR  THE  TECHSAT  21  MISSION 

Benjamin  J.  Bernocco 

Introduction 

The  Air  Force  is  currently  looking  for  new  technologies  to  take  it  into  the  twenty-first  century.  One  of 
those  technologies  high  on  the  Air  Force’s  list  is  Space-Based  Radar  (SBR).  SBR  will  provide  the  Air  Force  with  a 
virtual  presence  that  will  allow  constant  monitoring  of  situations  around  the  globe.  It  will  be  able  to  provide  theater¬ 
wide,  near  real-time,  undeniable,  all-weather,  and  round-the-clock  surveillance  of  any  place  on  the  globe.  The 
TechSat  21  program  was  setup  to  help  make  the  concept  of  SBR  a  reality  for  the  twenty-first  century.  It  will  consist 
of  clusters  of  small,  lightweight  formation  flying  satellites.  These  satellite  clusters  will  form  a  virtual  satellite  that 
functions  cooperatively  as  a  multiple  aperture  sparse  array.  This  array  would  be  capable  of  active  sensing, 
communications,  and  passive  radiometry.  Each  satellite  in  a  cluster  could  receive  its  own  as  well  as  other  satellites’ 
radar  pulses.  The  radar  pulses  are  of  different  frequencies  and  each  satellite  performs  range/Doppler  and  angle  of 
arrival  processing  on  the  pulses  [1].  In  order  to  maximize  data  collection  time  of  radar  signals,  the  thruster  firing 
time  of  the  satellite  must  be  minimized.  The  optimization  of  finite-thrust  spacecraft  trajectories  can  be 
accomplished  by  a  method  using  direct  collocation  with  nonlinear  programming. 

Methodology 

This  study  deals  with  the  optimization  of  finite-thrust  spacecraft  trajectories  by  using  the  direct  collocation 
with  nonlinear  programming  (DCNLP)  approach.  A  complete  history  of  the  development  of  the  DCNLP  approach 
can  be  found  in  references  2  and  3.  The  DCNLP  approach  approximates  the  spacecraft  trajectory  with  piecewise 
polynomials  that  are  represented  by  state  and  control  variables  at  a  given  number  of  nodes.  State  variables  are 
stored  in  the  state  vector  x  and  the  control  variables  are  stored  in  the  control  vector  u  .  The  spacecraft  trajectory  is 
divided  up  into  a  series  of  thrust  arcs  and  coast  arcs  and  the  thrust  arcs  are  subdivided  into  evenly  spaced  segments. 
Thrust  arcs  and  coast  arcs  require  different  means  of  evaluation.  In  a  thrust  arc  segment  each  state  variable  has  a 
state  trajectory  which  is  taken  to  be  a  Hermite  cubic.  This  is  a  “unique  cubic  that  goes  through  the  endpoints  of  the 
thrust  segment  with  the  appropriate  derivatives  that  are  dictated  by  the  evaluation  of  the  differential  equations  of 
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motion  at  the  endpoints”  [2],  Each  control  variable  is  assumed  to  have  a  linear  control  trajectory  across  the  thrust 
segment.  At  the  center  of  each  thrust  segment  a  collocation  point  is  found  and  then  the  difference  between  the 
derivative  of  the  Hermite  cubic  at  the  collocation  point  and  the  trajectory  equations  of  motion  evaluated  at  the 
collocation  point  is  calculated.  This  difference  is  called  the  defect  and  it  is  one  of  the  constraints  that  the  objective 
function,  the  thruster  firing  time,  is  subject  to.  When  the  defect  is  zero  the  differential  equations  of  motion  are 
satisfied  at  the  collocation  point  and  the  endpoints.  Allowing  the  trajectory  equations  of  motion  to  be  given  by 
x  =  f(x,u)  and  the  length  of  time  of  a  thrust  segment  to  be  T,  the  Hermite  interpolated  state  vector  at  the  collocation 
point  is  given  by 

xc  =  (1/2M*,  +  x,)  +  (T/8M/(*,,  m,)  -  f(x„  wr)]  (1) 

where  X\  and  xr  are  the  state  vectors  at  the  left  and  right  nodes.  The  linearly  interpolated  control  vector  is  given  by 

«c  =  (l/2)-(«,  +  ttr)  (2) 

where  «i  and  Mrare  the  control  vectors  at  the  left  and  right  nodes.  The  derivative  of  the  Hermite  interpolated  state 
vector  at  the  collocation  point  is 

xc  =  -[3/(2  T)].(X|  -  xT)  -  (1/4)-  [/(*,,  M|)  +  f{xt,  Mr)]  (3) 

Finally,  the  defect  vector  is 

d  =  f(xc,uc)  -  xc  (4) 

where  both  the  equations  of  motion  and  the  derivative  of  the  Hermite  interpolated  state  vector  are  evaluated  at  the 
collocation  point  [2], [3]. 

For  the  case  of  two-dimensional  motion,  which  this  study  focuses  on,  there  are  five  state  variables  and  one 
control  variable.  Based  on  this  information  there  will  be  five  equations  of  motion  for  the  system.  In  polar 
coordinates  using  canonical  units  (  p=  1.0)  the  two-dimensional  equations  of  motion  are 


*,  =  *3 

(5a) 

x2  =  X^l  X\ 

(5b) 

x3  =  *42/  *1  -  1  /  X\  +  *5  -sinCMj) 

(5c) 

x4  =  -XyX 4  /  X\  +  XS  -COS(Wi) 

(5d) 

x5  =  Xs2Ic 

(5e) 
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where  x\  is  the  orbital  radius,  x2  is  the  true  anomaly,  x3  is  the  velocity  component  along  the  radial  direction,  x 4  is  the 
velocity  component  along  the  local  horizontal  direction,  x$  is  the  magnitude  of  the  thrust  acceleration,  u\  is  the  thrust 
vector  angle  measured  from  the  local  horizontal  direction,  and  c  is  the  rocket  exhaust  velocity  [2], [4], [5]. 

A  coast  arc  is  handled  differently  from  a  thrust  arc  since  its  solution  can  be  determined  analytically.  For  a 
coast  arc,  integrals  of  the  motion  are  calculated  at  the  left  and  right  endpoints  of  the  arc  and  then  subtracted  to  yield 
a  generalized  defect.  When  the  generalized  defect  is  zero,  the  integrals  of  motion  at  the  left  and  right  endpoints  of 
the  coast  arc  are  equal.  For  the  case  of  two-dimensional  motion  four  integrals  of  the  motion  are  needed.  The 
integrals  of  motion  chosen  for  this  study  are  angular  momentum,  specific  energy,  longitude  of  periapsis,  and  the  fact 
that  the  thrust  acceleration  must  be  zero  across  the  coast  arc  [2].  The  two-dimensional  integrals  of  the  motion  in 


polar  coordinates  using  canonical  units  ( |i  =  1 .0)  are 

Qi(x)=x\-x4  (6a) 

Q2(x)  =  (x32  +  x?)  /  2.0  -  1.0  /  X\  (6b) 

Q3(at)  =  x2  +  tan'![  -xyxj  (  at42  —1.0/  *i)]  (6c) 

Q4(x)  =  x5  (6d) 


where  Qi(x)  is  the  angular  momentum  of  the  orbit,  Q2(x)  is  the  specific  energy  of  the  orbit,  Q3(x)  is  the  longitude  of 
periapsis  of  the  orbit,  and  Q4(at)  is  the  thrust  magnitude  [2], [4], [5].  The  generalized  defect  for  the  integrals  of  the 
motion  is  written  as 

q  =  Q(*i)-Q(*r)  (7) 

Besides  the  defect  and  generalized  defect  vectors,  equality  and  inequality  constraints  can  also  be  added  to  the 
analysis.  These  are  totally  up  the  user  and  can  include  things  like  constraints  on  the  initial  or  final  conditions  or  path 
constraints.  All  of  the  defects  and  constraints  are  used  as  bounds  on  the  objective  function,  what  you  want  to 
minimize.  In  the  case  of  spacecraft  trajectories  the  objective  function  is  the  sum  of  the  thrust  arc  durations.  The 
state,  control,  and  defect  variables  for  both  thrust  and  coast  arcs  plus  any  other  problem  constraints  make  up  a 
nonlinear  programming  problem.  To  solve  this  problem  the  variables  must  first  be  assembled  into  state  and 
constraint  vectors.  The  state  and  control  vectors  and  the  thrust  arc  duration  time  are  assembled  into  the  nonlinear 
programming  (  NLP)  state  vector  as  follows 

*T  =  [  x,T,  WlT,  *2t,  M2t,  . . *nT,  *nT,  tb  t2,  . . .,  tk]  (8) 
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where  X  is  the  NLP  state  vector,  n  is  the  number  of  nodes  that  the  trajectory  is  divided  into,  t  is  the  thrust  arc 
duration  time,  and  k  is  the  number  of  thrust  arcs  on  the  trajectory  [2].  The  defects,  generalized  defects,  equality 
constraints,  and  inequality  constraints  are  assembled  into  the  NLP  constraint  vector  as  follows 

CT=  [  db  q,,  d2,  q2,  dm,  qp,  h1(  h2, hj,  g1(  g2, .... gj  (9) 

where  C  is  the  NLP  constraint  vector,  m  is  the  number  thrust  arc  segments,  p  is  the  number  of  coast  arcs,  h  are  the 
equality  constraints,  j  is  the  number  of  equality  constraints,  g  are  inequality  constraints,  and  i  is  the  number  of 
inequality  constraints  [2].  The  NLP  state  vector  must  always  be  larger  in  length  than  the  NLP  constraint  vector  or 
an  optimized  solution  will  not  be  found.  Once  the  NLP  state  and  constraint  vectors  are  found,  a  solution  can  be 
computed  by  placing  these  vectors  into  an  optimization  subroutine. 

For  this  study  the  optimization  subroutine  used  was  the  Sequential  Quadratic  Programming  (SQP) 
subroutine  [6].  This  subroutine  approximates  the  nonlinear  problem  with  a  sequence  of  quadratic  programming 
problems.  It  minimizes  the  objective  function  subject  to  making  the  defects  and  constraints  values  as  small  as 
possible.  The  SQP  subroutine  must  be  provided  with  the  initial  state  and  control  conditions  and  with  subroutines 
that  evaluate  the  trajectory  equations  of  motion,  the  Hermite  cubic  and  its  derivative,  the  thrust  arc  defects,  the 
generalized  defects,  any  constraints  included,  and  the  objective  function.  It  also  requires  the  size  of  the  NLP  state 
and  constraint  vectors  as  a  check  to  make  sure  that  the  problem  is  solvable. 

The  two-dimensional  TechSat  21  example  was  setup  using  the  following  orbital  parameters  in  canonical 
units  (  p  =  1.0):  1)  the  orbital  period  equals  271 TU;  2)  the  semi-major  axis  equals  1.1098  DU/TU;  3)  the  eccentricity 
equals  0.0096;  4)  the  Earth’s  gravitational  constant  at  sea  level  equals  1.00  DU/TU2;  and  5)  the  rocket  engine 
specific  impulse  equals  0.37183  TU.  One  DU  equals  6378.1  km  and  one  TU  equals  806.81 1  seconds.  The 
trajectory  was  divided  into  twelve  thrust  arcs  and  twelve  coast  arcs  with  each  thrust  arc  divided  into  four  equa^ 
segments.  This  arrangement  occurs  over  a  two-year  period  and  results  in  a  trajectory  composed  of  sixty  nodes.  The 
equality  constraints  used  required  that  the  initial  and  final  eccentricity,  angular  momentum,  and  specific  energy 
calculated  from  the  state  and  control  variables  must  equal  the  actual  values  for  the  orbit  described  above.  No 
inequality  constraints  were  used  for  this  study.  The  initial  conditions  for  the  orbit  were  based  on  a  perfectly 
Keplerian  orbit  and  the  thrust  level  and  angle  were  randomly  guessed  at.  Results  based  on  this  setup  follow. 


1-6 


Results 


The  initial  conditions  in  Table  1  below  were  input  into  the  SQP  subroutine,  which  then  optimized  the  state 
and  control  values  to  yield  a  minimum  thruster  firing  time. 


Table  1:  Trajectory  Initial  Conditions  for  the  SQP  Subroutine 


Orbit 

Radius 

(DU) 

True 

Anomaly 

(rad.) 

Thrust 

(kg- 

DU/TU2) 

Thrust 

Angle 

(rad.) 

Thrust 

Firing 

Time 

<TU) 

Orbit 

Radius 

(DU) 

True 

Anomaly 

(rad.) 

Thrust 

(kg- 

DU/TU2) 

Thrust 

Angle 

(rad.) 

Thrust 

Firing 

Time 

(TU) 

1.088648 

0.0000 

8.00E-32 

-0.5236 

0.1 

1.088648 

6.2832 

8.00E-32 

-0.5236 

0.1 

1.088868 

0.2094 

8.00E-32 

-0.5236 

0.1 

1.088868 

6.4926 

8.00E-32 

-0.5236 

0.1 

1.089527 

0.4189 

8.00E-32 

-0.5236 

0.1 

1.089527 

6.7021 

8.00E-32 

-0.5236 

0.1 

1.090626 

0.6283 

8.00E-32 

-0.5236 

0.1 

1.090626 

6.9115 

8.00E-32 

-0.5236 

0.1 

1.092055 

0.8378 

0.00E+00 

0 

0 

1.092055 

7.1209 

0.00E+00 

0 

0 

1.093814 

1.0472 

8.00E-32 

-0.5236 

0.1 

1.093814 

7.3304 

8.00E-32 

-0.5236 

0.1 

1.095902 

1.2566 

8.00E-32 

-0.5236 

0.1 

1.095902 

7.5398 

8.00E-32 

-0.5236 

0.1 

1.097991 

1.4661 

8.00E-32 

-0.5236 

0.1 

1.097991 

7.7493 

8.00E-32 

-0.5236 

0.1 

1.100189 

1.6755 

8.00E-32 

-0.5236 

0.1 

1.100189 

7.9587 

8.00E-32 

-0.5236 

0.1 

1.102388 

1.8850 

0.00E+00 

0 

0 

1.102388 

8.1681 

0.00E+00 

0 

0 

1.104366 

2.0944 

8.00E-32 

-0.5236 

0.1 

1.104366 

8.3776 

8.00E-32 

-0.5236 

0.1 

1.106235 

2.3038 

8.00E-32 

-0.5236 

0.1 

1.106235 

8.5870 

8.00E-32 

-0.5236 

0.1 

1.107664 

2.5133 

8.00E-32 

-0.5236 

0.1 

1.107664 

8.7965 

8.00E-32 

-0.5236 

0.1 

1.108873 

2.7227 

8.00E-32 

-0.5236 

0.1 

1.108873 

9.0059 

8.00E-32 

-0.5236 

0.1 

1.109532 

2.9322 

O.OOE+OO 

0 

0 

1.109532 

9.2153 

0.00E+00 

0 

0 

1.109752 

3.1416 

8.00E-32 

-0.5236 

0.1 

1.109752 

9.4248 

8.00E-32 

-0.5236 

0.1 

1.109532 

3.3510 

8.00E-32 

-0.5236 

0.1 

1.109532 

9.6342 

8.00E-32 

-0.5236 

0.1 

1.108873 

3.5605 

8.00E-32 

-0.5236 

0.1 

1.108873 

9.8437 

8.00E-32 

-0.5236 

0.1 

1.107664 

3.7699 

8.00E-32 

-0.5236 

0.1 

1.107664 

10.0531 

8.00E-32 

-0.5236 

0.1 

1.106235 

3.9794 

0.00E+00 

0 

0 

1.106235 

10.2625 

0.00E+00 

0 

0 

1.104366 

4.1888 

8.00E-32 

-0.5236 

0.1 

1.104366 

10.4720 

8.00E-32 

-0.5236 

0.1 

1.102388 

4.3982 

8.00E-32 

-0.5236 

0.1 

1.102388 

10.6814 

8.00E-32 

-0.5236 

0.1 

1.100189 

4.6077 

8.00E-32 

-0.5236 

0.1 

1.100189 

10.8909 

8.00E-32 

-0.5236 

0.1 

1.097991 

4.8171 

8.00E-32 

-0.5236 

0.1 

1.097991 

11.1003 

8.00E-32 

-0.5236 

0.1 

1.095902 

5.0265 

0.00E+00 

0 

0 

1.095902 

11.3097 

0.00E+00 

0 

0 

1.093814 

5.2360 

8.00E-32 

-0.5236 

0.1 

1.093814 

11.5192 

8.00E-32 

-0.5236 

0.1 

1.092055 

5.4454 

8.00E-32 

-0.5236 

0.1 

1.092055 

11.7286 

8.00E-32 

-0.5236 

0.1 

1.090626 

5.6549 

8.00E-32 

-0.5236 

0.1 

1.090626 

11.9381 

8.00E-32 

-0.5236 

0.1 

1.089527 

5.8643 

8.00E-32 

-0.5236 

0.1 

1.089527 

12.1475 

8.00E-32 

-0.5236 

0.1 

1.088868 

6.0737 

0.00E+00 

0 

0 

1.088868 

12.3569 

0.00E+00 

0 

0 
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This  data  was  plugged  into  the  SQP  subroutine  along  with  the  subroutines  needed  for  the  DCNLP  calculations  and 
after  only  two  iterations  it  converged  to  an  optimized  solution.  The  optimized  trajectory  values  are  shown  in  Table 


Table  2:  Optimized  Trajectory  Values 


Orbit 

Radius 

(DU) 


True 

Anomaly 

(rad.) 


Thrust  Thrust 
(kg-  Angle 
DU/TU2)  (rad.) 


Orbit  True 

Radius  Anomaly 
(DU)  (rad.) 


Thrust  Thrust 
(kg-  Angle 
DU/TU2)  (rad.) 


1.088648 

0.0880 

2.11E-05 

0.4363 

0.1568 

1.096672 

6.4667 

2.47E-05 

0.4363 

0.1350 

1.089527 

0.2240 

2.11E-05 

0.4363 

0.1568 

1.096892 

6.5568 

2.49E-05 

0.4363 

0.1350 

1.090736 

0.3599 

2.11E-05 

0.4363 

0.1568 

1.097221 

6.6469 

2.50E-05 

0.4363 

0.1350 

1.092165 

0.4958 

2.09E-05 

0.4363 

0.1568 

1.097771 

6.7369 

2.47E-05 

0.4363 

0.1350 

1.093814 

0.6317 

0.00E+00 

0 

0 

1.098431 

6.8270 

0.00E+00 

0 

0 

1.104036 

1.3169 

2.11E-05 

0.4363 

0.1295 

1.104146 

7.5420 

2.48E-05 

0.4363 

0.1518 

1.105026 

1.3953 

2.06E-05 

0.4363 

0.1295 

1.105026 

7.6671 

2.47E-05 

0.4363 

0.1518 

1.106125 

1.4737 

2.19E-05 

0.4363 

0.1295 

1.106015 

7.7922 

2.47E-05 

0.4363 

0.1518 

1.107334 

1.5521 

2.33E-05 

0.4363 

0.1295 

1.107114 

7.9172 

2.45E-05 

0.4363 

0.1518 

1.108433 

1.6305 

0.00E+00 

0 

0 

1.108104 

8.0423 

0.00E+00 

0 

0 

1.107554 

2.1322 

2.28E-05 

0.4363 

0.1990 

1.112061 

8.2113 

2.47E-05 

0.4363 

0.2566 

1.110082 

2.3561 

2.31E-05 

0.4363 

0.1990 

1.114809 

8.5557 

2.47E-05 

0.4363 

0.2566 

1.11228 

2.5800 

2.29E-05 

0.4363 

0.1990 

1.117007 

8.9002 

2.46E-05 

0.4363 

0.2566 

1.114039 

2.8039 

2.21E-05 

0.4363 

0.1990 

1.118656 

9.2447 

2.46E-05 

0.4363 

0.2566 

1.115358 

3.0278 

0.00E+00 

0 

0 

1.119645 

9.5892 

0.00E+00 

0 

0 

1.11272 

3.2260 

2.30E-05 

0.4363 

0.1852 

1.119315 

9.3690 

2.21E-05 

0.4363 

0.2199 

1.11327 

3.4210 

2.3  IE-05 

0.4363 

0.1852 

1.118766 

9.6358 

2.94E-05 

0.4363 

0.2199 

1.11338 

3.6160 

2.34E-05 

0.4363 

0.1852 

1.117557 

9.9044 

1.92E-05 

0.4363 

0.2199 

1.11316 

3.8110 

2.30E-05 

0.4363 

0.1852 

1.115688 

10.1686 

1.64E-05 

0.4363 

0.2199 

1.1125 

4.0059 

0.00E+00 

0 

0 

1.11349 

10.4348 

0.00E+00 

0 

0 

1.111841 

4.1750 

2.46E-05 

0.4363 

0.1895 

1.109313 

10.4393 

0.00E+00 

0 

0.1834 

1.110632 

4.3789 

2.49E-05 

0.4363 

0.1895 

1.107224 

10.6302 

5.47E-06 

0.4363 

0.1834 

1.109313 

4.5827 

2.40E-05 

0.4363 

0.1895 

1.105026 

10.8213 

1.03E-06 

0.4363 

0.1834 

1.107994 

4.7866 

2.50E-05 

0.4363 

0.1895 

1.102717 

11.0124 

1.00E-06 

0.4363 

0.1834 

1.106675 

4.9904 

0.00E+00 

0 

0 

1.100409 

11.2035 

0.00E+00 

0 

0 

1.098101 

5.4931 

2.47E-05 

0.4363 

0.1114 

1.093484 

11.7332 

5.27E-07 

0.4363 

0.1162 

1.098431 

5.5336 

2.48E-05 

0.4363 

0.1114 

1.092495 

11.7843 

5.27E-07 

0.4363 

0.1162 

1.09887 

5.5740 

2.48E-05 

0.4363 

0.1114 

1.091506 

11.8354 

5.27E-07 

0.4363 

0.1162 

1.09942 

5.6145 

2.48E-05 

0.4363 

0.1114 

1.090626 

11.8867 

5.27E-07 

0.4363 

0.1162 

1.09986 

5.6549 

O.OOE+OO 

0 

0 

1.089857 

11.9377 

0.00E+00 

0 

0 

As  can  be  seen  from  Table  1  and  2,  the  optimized  trajectory  is  close  to  the  original  perfect  Keplerian  orbit.  This 
result  is  true  because  the  thrusters  are  only  firing  to  maintain  the  same  orbit;  to  cancel  out  the  general  perturbation 
included  in  the  equations  of  motion  of  the  system  via  the  thrust  acceleration  terms.  Figure  1  below  shows  a  plot  of 
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the  trajectory  initial  conditions  and  the  final  optimized  trajectory  values.  The  Earth  is  not  shown  to  scale  on  this 
figure. 


Figure  1:  Initial  and  Optimized  Trajectories 

The  major  differences  between  the  initial  and  the  optimized  values  are  the  thrust  magnitude,  the  thrust 
angle,  and  the  thruster  firing  time.  For  the  thrust  level  the  optimization  code  gave  values  which  stayed  between  0.1 
to  10  milli-Newtons  of  thrust.  At  one  point  in  the  trajectory  it  erased  a  thrust  segment  to  make  a  longer  coast  arc. 
The  thrust  angle  was  changed  to  a  fix  value  of  0.4363  radians  for  all  thrust  segments.  This  seems  odd  that  the  thrust 
segments  would  all  require  the  same  thrust  angle  but  no  error  was  found  in  the  coding  of  the  problem.  It  may  have 
something  to  do  with  the  way  the  optimization  subroutine  handles  the  control  values;  this  problem  will  have  to  be 
looked  at  in  greater  detail.  The  thruster  firing  time  increased  from  its  initial  value  but  this  was  necessary  in  order  for 
the  objective  function  to  meet  the  defect  and  equality  constraints  placed  on  it.  After  the  first  run  the  optimized 
trajectory  values  were  placed  back  into  the  SQP  subroutine  to  produce  new  optimized  values.  The  second  time 
through  the  routine  finished  after  only  one  iteration.  It  gave  the  exact  same  solution  as  is  shown  in  Table  2.  This 
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process  was  repeated  five  times  in  a  row  and  each  time  yielded  the  same  solution  as  in  Table  2  after  one  iteration. 
Such  consistency  in  getting  back  the  same  answer  implies  that  this  solution  is  a  fairly  optimized  solution  to  the 
problem  and  it  could  be  used  for  real  mission  planning.  Next  the  initial  conditions  were  offset  slightly  from  the 
original  values  and  it  converged  to  basically  the  same  solution  with  only  the  radius  and  the  true  anomaly  values 
changing  slightly.  A  really  bad  guess  for  the  trajectory  was  attempted  but  the  code  couldn’t  converge  to  a  solution. 
The  SQP  may  not  be  robust  enough  for  bad  guesses  so  a  more  efficient  and  robust  subroutine  like  NPSOL  or 
MINOS  should  be  used  if  the  there  is  no  good  guess  for  the  optimized  solution. 

Now  that  the  two-dimensional  code  is  working,  a  three-dimensional  code  needs  to  be  developed  since  all 
real  orbits  are  three-dimensional.  Currently  the  three-dimensional  equations  of  motion  and  integrals  of  the  motion 
are  developed  and  only  need  to  be  coded  into  the  DCNLP  program.  Once  the  three-dimensional  code  is  fully 
functional  perturbations  due  to  atmospheric  drag,  luni-solar,  and  the  Earth’s  oblateness  effect  can  be  included  as 
well  as  adding  multiple  satellites  to  make  a  realistic  model  of  the  TechSat  21  mission.  Numbers  returned  by  a  fully 
developed  three-dimensional  code  with  perturbations  could  be  used  for  actual  future  mission  planning. 

Conclusion 

An  optimal  finite-thrust  two-dimensional  spacecraft  trajectory  for  the  TechSat  21  mission  was  found  using 
the  direct  collocation  with  nonlinear  programming  approach.  This  approach  has  proven  itself  to  be  quite  effective  at 
calculating  optimized  trajectories.  The  optimized  trajectory  was  found  to  be  very  close  to  the  estimated  trajectory, 
which  instills  confidence  that  the  solution  is  a  valid  one.  Future  work  includes  expanding  the  code  to  three- 
dimensions,  the  inclusion  of  perturbations  caused  by  atmospheric  drag  and  the  Earth’s  oblateness  effects,  and  the 
inclusion  of  multiple  satellites  flying  in  formation.  The  only  improvement  to  be  made  is  to  get  a  more  robust  and 
efficient  optimization  subroutine  like  NPSOL  or  MINOS  that  can  handle  bad  guesses  for  the  initial  conditions.  This 
would  be  a  great  help  in  situations  where  no  obvious  solution  is  known  and  would  allow  greater  flexibility  in  the 
choice  of  initial  conditions. 
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Boulder,  CO  80309-0429 

ABSTRACT.  Perturbation  methods  for  linear  operators  are  commonly  used  in  the  analysis  of  systems  that 
tend  to  deviate  linearly  from  a  given  reference  model.  The  construction  of  an  operator  with  a  first  order 
perturbation  is  investigated  and  the  resulting  eigenvalue  series  is  constructed.  It  is  shown  that  a  simplified 
perturbation  series  can  be  obtained  for  matrix  operators  with  special  structure.  The  given  theory  is  applied 
to  a  reduced  order  model  (ROM)  control  scenario  and  an  algorithm  for  computing  an  0(e3)  eigenvalue 
approximation  is  described. 


Introduction 

Modem  modeling  of  large  dimensional  physical  systems  primarily  employs 
matrix  methods,  often  via  finite  elements  or  partial  differential  equations.  Due  to 
implementation  constraints  from  a  control  standpoint,  however,  an  adequate  model 
reduction  is  needed.  Model  approximations  can  be  used  to  design  controllers,  but 
compensation  for  observation  spillover  may  be  required  for  stability.  It  has  been  shown 
in  [1]  that  addition  of  a  residual  mode  filter  to  the  reduced  order  model  controller  can 
stabilize  the  entire  system  in  modal  coordinates.  The  question  now  becomes:  What 
components  of  the  spectrum  were  destabilized?  An  eigenvalue  perturbation  method, 
developed  in  [3],  has  been  proposed  for  full  rank  operators  with  no  multiplicities.  In  this 
paper,  we  have  expanded  this  technique  to  consider  multiple  eigenvalues  and  also  provide 
proof  of  the  order  of  approximation. 

Applied  mathematicians,  engineers  and  physicists  frequently  use  perturbation 
methods  to  solve  a  variety  of  different  problems.  Even  though  solutions  to  intractable 
nonlinear  differential  equations  can  be  approximated  using  perturbation  theory,  this 
paper’s  focus  will  only  encompass  linear  deviations  of  finite-dimensional  linear 
operators.  The  reader  is  also  referred  to  an  “in-depth”  treatment  of  the  perturbation  theory 
of  linear  operators  [4],  where  the  subsequent  derivations  become  special  cases. 

The  following  definitions  lay  the  groundwork  for  the  construction  of  an 
eigenvalue  perturbation  series.  The  type  of  structure  imparted  upon  a  matrix  operator  of 
interest  is  illustrated  and  the  corresponding  eigenvalue  series  is  given.  An  example  of  the 
application  to  the  control  design  of  a  spectral  system  is  explained.  A  proof  and  the 
numerical  algorithm  for  this  approach  are  given  in  the  appendices. 

Operator  Perturbations 
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(1) 


Beginning  with  an  operator  with  a  first  order  perturbation, 

A(£)  =  A°+£AA 

A0  e  RLxL  (diagonalizable,  full  rank) 

£ e  C,  M  e  RLxL 

the  resolvent  of  A(s)  is  defined  by  the  following  inverse: 

R(X,A{e))  =  [A{£)-XLY 
and  is  analytic  everywhere  except  at  the  eigenvalues  of  A(£).  It  has  been  shown  in  [3] 
that  the  resolvent  can  be  written  in  series  form, 


R(A,  A(£))  =  R(X,  A0 )  +  X*(n)  (X)£n 

n  *  1 

R{n)  (A)  m  R(A,  A°)[-AA  R(A,  A0 )]" 
and  the  series  is  convergent  if  the  following  inequality  holds. 


o<||a^(/M°)| 

The  spectral  factorization  of  A0  can  be  expressed 
projections,  Pt. 


J|f|<  1 

as  a  sum  of  S  orthogonal  eigen- 


(3) 


A°  =UAU~'  =  iv,'  =  t*,VE,U-'  (4) 

*=1  /=1  /-I 

'0  0  .  0" 
o  I„x.  •  0 

0  0  .0 

Note  that  S  <  N  if  any  of  the  eigenvalues  in  A  are  repeated.  The  multiplicity  of  the  /  * 
eigenvalue  is  computed  by  taking  the  trace  of  the  /  *  factor  in  equation  (4). 


Tr[u  E,U~l\=  Tr[E,  lT'u\=  Tr[E,  ]  =  m, 

Similarly,  the  resolvent  can  be  written  as  a  sum  of  projections: 

m,A°)='ZT^T  (5) 

i~\  Aj  —  A. 

and  the  reduced  resolvent,  Sn(/l),  is  defined  in  equation  (6)  as  the  complement  of  the  n  * 
projection. 

(«> 

l-l  —  A 
l^n 

Each  of  these  eigen-projections  can  be  recovered  from  the  integration  of  the  resolvent 
around  a  closed  curve  in  the  complex  plane. 


Operator  Partitioning 


Assume  that  the  linear  operator  A(s)  can  be  partitioned  in  the  following  manner: 


A(s)  = 


0  A 2 


0  A41 


=  A0  +  s  A  A 
A0  e  XBD  (Block  Diagonal) 

A A  e  XqD  (Block  Off-Diagonal) 

The  first  factor,  A0,  is  an  element  in  the  space  of  block  diagonal  matrices.  The  perturbing 
matrix,  AA,  is  a  member  of  the  complementary  block  off-diagonal  space.  In  fact,  matrix 
multiplication  of  these  operators  can  be  interpreted  as  mappings  between  these  spaces. 

It  is  said  that  A  and  B  are  partitioned  when  A  is  a.  complement  matrix  of  B  (A  and 
B  do  not  overlap)  and  A  e  XBD;  B  e  XqD,  then  the  following  definitions  hold: 

(1.)  A"  s  XBD  V  n  =  1,2, ... 


(2  .)Bne 


XBD  for  n  even 
X0D  f°r  n  °dd 


A-B]  „ 
(3.)  \<=  Xc 

v  JB-A\  c 


Eigenvalue  Perturbations 

The  if1  projection  of  an  operator,  A{s),  can  be  written  as  the  corresponding 
eigenvalue  scaling  that  projection. 

A(S)  Pk  (£)  =  Pk  (£)  A(e)  =  ;.k  {£)  Pk  (£)  (9} 

Pfe)  is  determined  from  integrating  the  resolvent,  R(Z,  A(e)),  around  the  kP  eigenvalue 
and  results  in  an  infinite  series  in  s. 

P(£)  =  — =  +  ft *"”(•*)  « 

2  n  j  2  7t  j  r,  "=i 

=  Pk  +  X7*"’ £*  (10) 

n- 1 

By  taking  the  trace  of  equation  (9),  the  A*  eigenvalue  series  is  recovered: 

Zk(e)=—Tr[A(e)Pk{e)}  (11) 

=  A  + 
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where  A0  represents  an  eigenvalue  of  A0.  If  A(s)  is  partitioned  into  the  form  given  in  (8), 
the  previous  equation  is  simplified  to  an  even  series. 

K  (*)  «  A  +  4"  e1  +  4"  s‘  + ...  =  4  +  (12) 

n=l 

Proof  of  this  result  is  given  in  Appendix  A. 

Determination  of  Instability  in  Linearly  Perturbed  Systems 

At  this  point,  the  application  of  this  method  to  the  state  space  control  approach 
given  in  [1],  [2],  and  [5]  is  considered.  Given  a  linear  plant  governed  by  the  following 
vector  differential  equation: 


x  =  Ax  +Bu 

~  ~  (13) 

y = Cx+Du 

x  =  plant  state  vector 
u  =  control  input  vector 
y  =  output  observation  vector 

the  objective  is  to  stabilize  (13)  through  the  use  of  M  actuator  inputs  and  P  sensor 
outputs.  Matrix  operators  B  and  C  will  have  the  finite  rank  of  M  and  P,  respectively,  with 
the  linear  differential  operator.  A,  possibly  of  infinite  dimension. 

Using  the  ROM  (reduced  order  model)  methodology,  it  is  possible  to  transform 
the  plant  to  a  modal  coordinate  system  and  separate  (13)  into  two  state  equations. 
Assuming  the  feed-through  term,  D,  is  zero: 


—  An  xN  +  B N  u_ 

yN  =  CN 


(14a) 


x_R  =ARx^R  +  BRu_ 
}Lr=CrXr 
y  -  y  +  y 


(14b) 


The  preliminary  control  design  of  the  entire  plant  is  based  upon  the  dynamics  of  (14a).  >  A 
finite  number  of  modes  are  chosen  for  AN,  which  is  an  operator  of  rank  N.  Equation 
(14b)  consists  of  modes  that  are  not  of  immediate  consideration  and  are  open  loop  stable. 
(Re  (A  (AjJ)  <  0)  The  use  of  state  estimation  and  N  state  feedback  results  in  the  system  of 
equations  in  (15). 


^ ,v  —  xn  +  Bn  u  +  Kn (y  y^) 

L = c*  ** 


(15) 


u  =  Gnxn 

Defining  the  error  between  estimated  and  actual  states  by  a  vector  eN,  the  differential 
equation  governing  error  progression  becomes  the  difference  between  state  equations: 
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(16) 


§.N  ~  In  —S 

eN  =  (  An  —  KnCn)  eN  +  KnCrxr 
If  controllability  and  observability  requirements  are  met  for  the  system,  arbitrary  decay 
rates  for  the  designed  subsystem  can  be  achieved  by  choosing  the  appropriate  estimator 
gain,  K„,  and  feedback  gain,  GN,  matrices  for  the  desired  eigenvalues  of  A-KC  and 

A+BG.  iii  c 

Writing  equations  (14a),  (14b),  and  (16)  in  matrix  form,  the  probable  causes  of 

instability  in  the  closed  loop  system  can  be  clearly  determined. 


• 

Is 

An+BnGn  BnGn  0 

*iV 

Iv 

- 

0  An  -KnCn  k„cr 

* 

£iV 

In 

B,Gb  Ar 

1 - 

\H 

>3 

Cross-coupling  terms  in  (17),  illustrated  in  bold,  characterize  the  introduction  of 
estimator  feedback  which  can  destabilize  the  residual  state  equation.  One  option  to 
remedy  this  problem  is  given  in  [1],  in  which  a  residual  mode  filter  (RMF)  is  developed 
to  compensate  for  eigenvalues  that  are  drawn  unstable.  However,  knowledge  of  those 
open  loop  eigenvalues  contained  in  AR  is  needed  in  order  to  implement  the  method. 
Previous  work  by  Gooyabadi  [3]  offers  a  limited  perturbation  solution  to  this  problem. 

Partitioning  equation  (17)  according  to  the  method  given  in  section  (2),  die  matrix 
Ac  can  be  considered  to  be  the  sum  of  a  stable,  linear  operator,  A0,  and  a  perturbing 
operator,  A A. 


A(s)  = 


an  +bngn 

BnGn 

0  ' 

0 

0 

0 

0 

”  KnCm 

0 

+  £ 

0 

0 

KnCr 

0 

0 

Ar_ 

_BrGn 

BrGn 

0 

(18) 


=  A°  +  £  A A 

Our  objective  is  to  find  out  which  eigenvalues  of  A0  are  driven  unstable  by  AA  through 
the  computation  of  the  first  few  non-zero  terms  of  the  perturbation  series  (12).  Taking  e 
=  1  would  give  equation  (17),  but  convergence  of  the  eigenvalue  series  solely  rests  upon 
the  condition  that  \\AA  R{X,  i4°)||  <  1.  Therefore,  it  is  possible  that  (12)  is  divergent  for 

arbitrary  \\A4  R(X,  A°)\\  *  1  and  may  have  asymptotic  convergence.  .  ■ 

An  example  of  implementation  of  this  perturbation  scheme  is  given  for  a  simple 

Euler-Bemoulli  beam  with  pinned  ends: 

mzSxA)-2Zzlxx{xA)  +  E-Iz^(xA)  =  u{t)8{x-QAA)  (19) 

z(0,t)=z(l,  r)  =  0 
B*C”  z  xx(0,  t)  =  z  „(/,  /)  =0 

y(t)=  fz(x,t)5(x-0.9-l)  dx  =  z(0.9-l,t) 

All  of  the  constants  in  (19)  are  normalized  with  the  exception  of  the  damping  ratio.  ( \Z\ 
5.0  1CT6 )  Actuation  of  the  beam  occurs  at  the  point  x  =  0.1.  Observations  of  the  beam 
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displacement,  y  (t),  are  taken  at  x  =  0.9.  Using  the  method  of  eigenfunction  expansions,  a 
temporal  differential  equation  can  be  constructed: 

zk(t)  =  -(k7ty  zk (t)-2%(k 7tf  z4(O  +  &(0.1)  u(t)  (20) 

The  spatial  eigenfunctions,  <f>k,  transform  the  input,  u  (t),  into  the  modal  system  of 
equations  given  in  (20).  Control  is  introduced  to  dampen  the  first  few  modes,  where  k  = 
1, 2,  and  3.  There  exist  an  infinite  number  of  residual  modes  ( k  >  4),  but  for 
computational  purposes  only  the  next  12  are  considered.  Putting  the  equations  derived 
from  (13)  in  the  form  given  in  (18),  we  can  calculate  an  estimate  of  the  residual 
eigenvalue  perturbation  resulting  from  the  control  design: 


Ms)  (2d 

The  results  of  this  O  {£)  approximation  are  illustrated  below  in  Figure  1 .  from 
this  graph,  it  can  be  determined  that  modes  n  =  5,  7  are  driven  unstable.  Using  these 
modes  in  an  RMF  control  design,  exponential  stability  for  the  entire  system  can  be 
achieved. 
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Conclusions 

A  generalized  approach  to  the  perturbation  method  outlined  in  [3]  has  been 
developed.  Given  a  block  diagonal  linear  operator  with  a  first  order  perturbation  in 
partitioned  form,  it  has  been  shown  that  there  exists  a  corresponding  even  perturbation 
series.  An  algorithm  is  outlined  to  find  the  O  {£)  correction  for  block  diagonal  operators 
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containing  eigenvalue  multiplicities.  Application  of  this  technique  to  a  simple  structural 
model  illustrates  that  this  method  is  a  viable  approach  to  ROM  control  design.  This 
procedure  is  shown  to  be  successful  in  other  structural  problems  with  no  repeated 
eigenvalues  using  finite  element  models  [3]. 
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Appendix  A. 


Equation  (12)  is  derived  from  (1 1)  as  follows: 

Let  A(s)  =  A°  +  sAA,  where  A0  e  Xbd  and  AA  e  Xod- 

Xk(e)=—Ti{A{e)Pk(e)\ 

mt 

= — T\{A{s)-??k\)  Pk  (£)+X°kPk  (£)} 

mk 

=  ??k+—T\{A{e)-A\\)Pk{e)\ 

mk 

\(s)-X°k  =-T?{(A(£)-All)Pk(£)} 

mk 

Rewriting  part  of  the  previous  result: 

(A(s ) -X°kI)Pk  (s)  =  ^~<f  (A(s) - X°k I )  R(X,  A{e))  dX 
2*J  t\ 

=  —  ^{A{£)-Xl  +  Xl-Xl\)R{X,A{£))dX 
2  *j  r* 

=  —  <f  ( A(£ ) -XI) R(X,  A(£))  dX  +  d*(A - X°k )  R(X,  A(£))  dX 

2jz  ]  |_r*  r. 

The  first  integral  in  (B2)  is  analytic  since  (A(£)  -  X  I)  R(X,  A(£))  =  I.  Therefore: 


X(£)-X°k  = 


—  Tr  <S(X-X \)  R{X,A{£))  dX 
J  m  k  Lr, 


A  series  formis  desired  for  (B3).  Rewriting  the  kernel  of  this  equation: 

(X-X°k)  R(X,  A(£))  =  (X-X°k)  R(X,  A°)  +  £(X-X°k)  R(n}  (X)  £n 

n=  1 

Integrating  the  first  term  in  the  right  hand  side  of  the  equation,  we  find  that  this  is 
completely  analytic  inside  Tk- 

-  A°t  )  R  (A,A°  )  dX  =  ja  -  A  )£-A—  M 

r  r\  I  -  A 


r*  * 


From  equation  (B4): 


X(s)-X°k= - - - Tr  cfY  {X- X°k)R("\X)£n  dX 

In  jmk  £ 


In  jmk 


■Tr  fz<-  1)"(/1-A®)/?(A,^0)[a^/?(A,^0)]”  £n  dX 
_r» 


Using  the  cyclic  permutation  property  of  the  trace  operator: 

Tr\—[AA  R(A,  A0  )]"1  =  n •  t\r{A,  A0 )  [a^  R(A,  A0 )]"  ]  (B6) 

dX  J 

Substitution  of  (B6)  into  equation  (B5)  yields  a  form  that  can  be  integrated  by  parts: 

In  jmk  |_r'  "  dA  J 

2njmk  ^  dA\JzA  n  J  r,  *»'  n 

The  first  term  in  (B7)  can  be  evaluated  about  the  circle  centered  at  A0.  Making  the 
substitution:  A  =  A\  +  r  e*  *,  this  integral  is  zero: 

2x 

S_±\yt}L(A-Al)[MR(A,A0)]£n  dA  =  ~ ~ [M R&  +  r e* > ' A° )1" £‘ "  =0 

fdA\jz t  J  »='  ”  o 

The  remaining  non-zero  term  in  (B7)  gives  a  series  in  s.  However,  the  goal  is  to  prove  that 
the  odd  terms  are  zero.  By  equation  (10),  it  is  seen  that  the  remainder  of  (B7)  can  be 
rewritten: 

"  n—\ 

Ak(s)-A°k  = - - — Tr  <S £(-A4 R(A, A0 ))•  ^—  [a4 R(A, A0 )]  s"  dA 

2njmk  [r*  « 

(  Y 

=  J-rr  AA-f—  —  cf(-l)”-1  R(A,  A0)[AA  RjAM^  dA 

™k  "='  n  2  ^  7  r,  R^\A)  '  J_ 


Tr  A4-£- 


oc  *  p(”~l) 


For  each  coefficient,  A(n\  in  this  expansion: 

=-l—Tr[AAP{kn-')]  (B8) 

mk  n 

Integrating  P  <n"1>  yields  an  expansion  in  terms  of  the  reduced  resolvent,  S  (kK),  and  the 
projection,  P  °. 

£  A  A  Slfl)(At)AA  S(k'>\Ak)...  A  A  Sj'-’OM 

/,  +  ...  +  /„  -  n  ~\ 

where  5 (0)  a  P  0  and  S (i)  =  5  '*  for  i  >  1 . 

Writing  (B9)  in  terms  of  an  ordered  multiplication: 
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4n)  =— Tr 
mkn 


I 

/,+...+/,=»- 1 
/^O 


fiM'-U)] 


(BIO) 


Since  A/4  e  Xod  and  S®  e  Xbd  for  ally  >  0  by  definition  (1)  of  section  (2),  then  A A 
S  0)  e  Xod  by  definition  (3).  Therefore,  equation  (BIO)  is  the  sum  of  a  multiplication  of  n 
block  off-diagonal  matrices  that  are  contained  in  either  Xod  or  Xbd,  depending  on  whether 
n  is  odd  or  even.  The  trace  of  a  block  off  diagonal  matrix  is  zero,  therefore  the  only  non¬ 
zero  coefficients  in  (1 1)  are  even.  Q.E.D. 
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Appendix  B. 


Assume  all  of  the  eigenvalues  of  the  full-rank  block  diagonal  matrix,  A0  e  RL  x  , 
and  the  corresponding  right  and  left  eigenvector  matrices,  U and  UA  can  be  found.  An 
algorithm  for  computing  the  first  two  non-zero  terms  in  the  eigenvalue  perturbation  series  is 
given  below. 

1.)  Diagonalize^0,  where  the  eigenvalues  in  A1’1  are  ordered  Re(  AMi)  >  Re(  X''i)  >  ...  > 
Re(  1UL): 


[>'  0  1 

‘l/u  0  " 

r-u 

A  0 

1 

o 

1 

M 

o 

_ 1 

_  0  U2'\ 

— 2,2 

_  °  A  J 

_  o  (u»y\ 

2.)  Form  the  eigen-projections P  yk  from P k  for  each  distinct  eigenvalue: 


A°  = 


:  1.1 


N\  , 

I A 


*= 1 


0 


N2-7  7 

I  tip* 


A- 1 


N(l)  =  #  of  eigenvalues  in  A 
N(2)  =  #  of  e.v.  in  A2,2. 

For  i  =  1  to  2, 
k=  1 
L  =  N(p 

P’  =  0  LxL 

R^CoMf/^Row^t/Uy1)  ** 
m(k,  i)  =  1 
For  j  =  2  to  L, 

Ptemp  =  Col^O  ®  R0Wj((f/ 

If  AMj.i  =  Awj, 

P’  =  P’  + Ptemp 
m(k,  i)  =  m(k,  i )  +  1 


Else, 


P'\  =  P’ 

X  l\  =  X  *’ j.i 
&  =  1 
m(£,  /)  =  1 
P  —  P temp 


End  If 

End 

P  Uk  =  P  ’ 

Auk=  Aul 

s<0  =  * 

End 

**  <g>  denotes  an  outer-product. 
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3.)  Compute  reduced  resolvents,  correcting  coefficients,  and  eigenvalue  estimates: 


For  k=  1  to  S(.l), 


5(2)  p2-2 

C2-2  -  V _ l _ 

*  Z- 1  il,1  i2,2 

7=1  Ak  A) 

UT  =-7rr:7>  I^1'2 

m(k9 1) 


For  k—  1  to  S(2), 


so)  Pu 

cw  -  Y  l _ 

*  La  i2,2  jl.l 

7=1  AJ 

te!f  =— !— 7>-|/Uu$,uA4ui,P] 

'  *  '  m(k,  2)  1 

32,2.  x_  i2,2  i  f j2,2 V2)  2 


A*(.£)  =  Ak  +1^  7  * 


End  If 
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PROFILE  IN  A  THIN  COMPOSITE  PART 


Jeff  M  Ganley 
Graduate  Research  Assistant 
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University  of  New  Mexico 


Abstract 


Composite  materials  are  advantageous  for  aerospace  applications  due  to  their  light  weight  and  high 
stiffhess/weight  ratio.  However,  composite  parts  which  are  autoclave  cured  experience  problems  associated 
with  cure  residual  stresses,  including  a  reduction  in  load  carrying  capacity  and  ‘spring-in’,  a  permanent 
deformation  due  to  the  residual  stresses  (Stover,  1993).  This  paper  develops  an  experimental  procedure  that  can 
be  used  to  obtain  the  residual  stress  profile  in  a  thin  composite  part.  This  procedure  involves  progressive  cutting 
of  the  composite  part  (i.e.  stress  relief)  and  associated  strain  gage  monitoring.  The  strain  gage  data  is  then 
combined  with  a  finite  element  analysis  to  determine  the  through-thickness  residual  stress  profile  in  the 
composite  part.  This  analysis  procedure  can  be  used  as  a  tool  in  the  understanding  and  prediction  of  spring-in, 
or  simply  for  measuring  the  residual  stresses  in  a  composite  part. 
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DETERMINATION  OF  THE  RESIDUAL  STRESS 
PROFILE  IN  A  THIN  COMPOSITE  PART 


Jeff  M  Ganley 


Introduction 

Composite  materials  are  advantageous  for  aerospace  applications  due  to  their  light  weight  and  high 
stiffness/weight  ratio.  Many  other  factors  give  composites  superior  performance  over  traditional  materials, 
including  the  fact  that  composite  materials  can  be  tailored  to  meet  specific  design  and  loading  requirements. 
However,  a  significant  problem  in  the  composites  industry  is  that  autoclave  cured  composite  parts  experience 
problems  associated  with  curing  residual  stresses,  including  a  reduction  in  load  carrying  capacity  and  spring-in. 
This  paper  will  present  an  experimental  methodology  for  determining  the  residual  stress  profile  in  a  thin 
composite  part.  This  methodology  was  originally  developed  as  a  tool  to  aid  in  understanding  spring-in  in  thin 
composite  parts.  However,  the  methodology  can  easily  be  extended  to  either  thick  composite  parts  or  isotropic 
material  applications.  In  fact,  the  methodology  is  applicable,  with  certain  limitations,  to  any  solid  where  the 
residual  stress  profile  is  desired. 

Determination  of  Residual  Stress  Profile 

The  purpose  of  the  methodology  developed  in  this  paper  is  to  be  able  to  precisely  determine  the  stress 
(or  strain,  given  a  linear  elastic  assumption)  profile  through  the  thickness  of  a  thin  composite  part.  For  most 
composite  materials,  including  the  IM7/977-2  carbon  fiber/epoxy  material  used  in  the  current  research,  the 
material  response  is  linear  elastic  to  initial  failure.  Thus,  the  linear  elastic  assumption  is  typically  a  good  one  for 
stresses  below  initial  failure  thresholds.  The  residual  stress  profile  determination  will  be  accomplished  through  a 
methodology  that  involves  a  step-by-step  cutting  of  the  composite  part  through  the  thickness,  combined  with 
finite  element  modeling  using  IDEAS  Master  Series  4.0  finite  element  software  (SDRC,  1990).  A  brief 
overview  of  the  procedure  for  finding  the  residual  stress  profile  through  the  composite  part  thickness  is  given  by 
the  following: 
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1.  Mount  strain  gages  to  the  composite  part  on  the  inside  and  outside  faces. 

2.  Cut  the  part  through  the  thickness  near  the  strain  gage  in  ten  equal  steps,  stopping  after  each  cut  to  read 
the  strain  gage  output.  The  strain  change  measured  by  the  gage  at  each  step  is  due  to  the  relieving  of 
residual  stresses  at  the  cut. 

3.  Model  the  part  in  IDEAS,  creating  finite  element  models  for  each  of  the  ten  cuts.  Place  a  unit  load  at 
each  modeled  cut  and  record  the  associated  strains  from  IDEAS  at  the  real  strain  gage  location. 

4.  Correlate  the  IDEAS  strains  with  the  actual  strain  gage  strains  to  obtain  the  residual  stresses  in  the  part  at 
each  cut. 

20  Laver  Part  Residual  Stress  Profile  Experiment 

The  procedure  given  above  was  completed  for  a  20  layer  thick,  17.8  N  winding  tension  IM7/977-2 
unidirectional  hoop  composite  part.  The  part  had  the  following  dimensions:  thickness  (radially)  =  0.3175  cm, 
width  =  4.064  cm,  and  original  inner  radius  =  6.585  cm.  In  preparation  for  the  test,  the  part  was  cut 
longitudinally  and  allowed  to  spring-in,  then  the  overlap  due  to  the  spring-in  was  cut  out.  The  part  was  then 
mounted  with  8  strain  gages,  mounted  2.54  cm  apart  center  to  center,  on  alternating  faces  of  the  part  (4  on  the 
inside  face  and  4  on  the  outside  face).  The  part  was  then  mounted  in  a  Bridgeport  C&C  machine  and  secured 
(i.e.  fixed  condition)  at  2.54  cm  below  the  first  strain  gage.  The  machine  was  then  programmed  to  cut  the  part  in 
ten  equal  steps  at  the  centerline  of  the  first  strain  gage  from  inside  to  outside  (the  first  gage  was  mounted  on  the 
outside  face).  After  each  of  the  ten  cuts,  the  machine  was  paused  for  approximately  30  seconds  while  the  output 
from  the  strain  gage  signal  conditioner  was  recorded.  When  the  cutting  process  at  the  first  gage  was  completed, 
the  part  was  rotated  to  align  the  cutting  blade  with  the  centerline  of  the  second  gage  and  re-clamped,  and  the 
cutting  process  repeated.  For  each  strain  gage,  the  cutting  was  from  the  opposite  face  of  the  part,  toward  the 
gage,  at  the  centerline  of  the  gage.  In  addition,  for  each  cutting  process  there  was  2.54  cm  of  free  composite  part 
above  the  cut,  and  2.54  cm  from  the  cut  location  to  the  fixed  condition.  This  process  was  repeated  for  all  eight 
strain  gages. 

The  reason  for  cutting  at  the  centerline  of  the  strain  gage,  which  will  cut  through  and  destroy  the  gage 
on  the  final  cut,  is  that  this  will  give  more  robust  results  than  mounting  the  strain  gage  offset  to  the  cut.  As  can 
be  seen  in  Figure  1  on  the  following  page,  the  FEM  strain  values  due  to  the  cutting  and  relieving  of  residual 
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stresses  in  the  part  are  fairly  constant  over  a  range  approximately  equal  to  the  thickness  of  the  part  (-5  to  5  in 
Figure  1).  However,  the  strains  quickly  drop  off  as  you  move  away  from  the  center  of  the  cut. 


Mat'l  x-strain  (Outside) 


Figure  1  -  Strains  on  the  Outer  Face  of  the  20  Layer  Part  Due  to  the  First  Cut 
as  a  Function  of  Distance  from  the  Cut 

If  the  strain  gage  is  placed  offset  to  the  cut  and  the  strain  gage  is  not  precisely  located  relative  to  the 
cut,  the  actual  gage  location  and  the  FEM  modeled  gage  location  will  not  be  the  same,  and  thus  the  two  will 
measure  differing  strains  (refer  to  Figure  1  above).  This  will  lead  to  error  in  the  prediction  of  the  residual 
stresses  in  the  part.  However,  if  the  strain  gage  is  placed  at  the  centerline  of  the  cut,  a  slight  misplacement  of  the 
gage  will  not  effect  the  results  significantly.  This  is  because  the  strains  near  the  cut  centerline  are  relatively 
constant  and  because  of  the  symmetry  of  the  strain  values  about  the  centerline  of  the  cut  (refer  to  Figure  1). 
Thus,  a  slight  misalignment  of  the  strain  gage  relative  to  the  cut  will  not  significantly  change  the  average  strain 
over  the  length  of  the  strain  gage. 

In  addition  to  the  requirements  given  above,  it  is  important  to  consider  other  factors  when  performing  a 
residual  stress  profile  test  on  a  thin  composite  part.  First,  one  must  be  sure  to  use  a  sealed  strain  gage.  If  not,  the 
carbon  dust  that  results  from  the  cutting  will  settle  on  the  face  of  the  strain  gage  and  short  the  gage,  yielding 
inaccurate  strain  readings.  Second,  because  of  the  thin  nature  of  the  part  and  the  placement  of  the  strain  gage 
close  to  the  cut,  the  strain  gage  and  part  will  experience  heating  due  to  the  cutting  operation.  This  must  be 
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accounted  for  during  the  testing.  This  was  accomplished  in  this  research  by  selecting  a  thin  diamond  cutting 
blade  (i.e.  less  heat  production),  and  allowing  the  part  and  gage  to  return  to  room  temperature  prior  to  recording 
the  strain  gage  output. 

The  results  of  the  20  layer  part  residual  stress  profile  test  are  given  in  Table  1  below.  The  values  in  the 
first  part  of  the  table  are  from  the  strain  gage  signal  conditioner  and  are  given  in  Volts.  The  second  part  of  Table 
1  gives  the  average  values  from  the  eight  strain  gages  in  the  first  part  of  the  table  (in  Volts,  V0),  then  converts 
the  average  to  microstrains  (pe)  using  the  following  equation: 

pe  =  V0  *  382.7751 196  (Equation  1) 


Table  1  -  Residual  Stress  Profile  Strain  Gage  Results  (20  Layer  Part) 


Cut 

1st  gage/ 
inside-out 

2nd  gage  / 
outside-in 

3rd  gage  / 
inside-out 

4th  gage  / 
outside-in 

5th  gage  / 
inside-out 

6th  gage  / 
outside-in 

7th  gage  / 
inside-out 

8th  gage  / 
outside-in 

No  Cut 

0.002 

0.004 

0.002 

0.000 

0.002 

0.000 

0.003 

0.000 

0.03175  cm 

-0.403 

-0.237 

-0.393 

-0.217 

-0.436 

-0.208 

-0.504 

-0.293 

0.06350  cm 

-0.942 

-0.600 

-0.952 

-0.614 

-0.942 

-0.578 

-1.060 

-0.713 

0.09525  cm 

-1.416 

-1.066 

-1.432 

-1.093 

-1.335 

-1.058 

-1.520 

-1.235 

0.12700  cm 

-1.660 

-1.620 

-1.682 

-1.669 

-1.521 

-1.626 

-1.756 

-1.835 

0.15875  cm 

-1.746 

-2.207  i 

-1.771 

-2.270 

-1.577 

-2.232 

-1.816 

-2.444 

0.19050  cm 

-1.727 

-2.743 

-1.750  : 

-2.825 

-1.546 

-2.813 

-1.770 

-2.970 

0.22225  cm 

-1.680 

-3.067  ! 

-1.716  i 

-3.160 

-1.520 

-3.200 

-1.690 

-3.280 

0.25400  cm 

-1.667 

-3.150  j 

-1.730 

-3.270 

-1.558 

-3.340 

-1.616 

-3.420 

0.28575  cm 

-1.916 

-3.120 

-1.947  | 

-3.420 

-2.055 

-3.300 

-1.819 

-3.630 

0.31750  cm 

2.025 

1.080  1 

-0.098 

1.411 

12.520 

-3.930 

-2.507 

1.500 

Cut 

Average  of 
1,3,5  &  7 

Average  of 
2,4,6  &  8 

Avg.  strain 

1.3, 5, 7  (pc) 

Avg.  strain 

2, 4, 6, 8  (pc) 

No  Cut 

0.0023 

0.0010 

0 

0 

0.03175  cm 

-0.4340 

-0.2388 

-166.12 

-91.39 

0.06350  cm 

-0.9740 

-0.6263 

-372.82 

-239.71 

0.09525  cm 

-1.4258 

-1.1130 

-545.74 

-426.03 

0.12700  cm 

-1.6548 

-1.6875 

-633.40 

-645.93 

0.15875  cm 

-1.7275 

-2.2883 

-661.24 

-875.89 

0.19050  cm 

-1.6983 

-2.8378 

-650.05 

-1086.22 

0.22225  cm 

-1.6515 

-3.1768 

-632.15 

-1215.98 

0.25400  cm 

-1.6428 

-3.2950 

-628.80 

-1261.24 

0.28575  cm 

-1.9343 

-3.3675 

-740.38  | 

-1289.00 

0.31750  cm 

- 

- 

! 

- 
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The  values  for  the  average  strain  given  in  Table  1  on  the  previous  page  will  be  correlated  with  the 
IDEAS  FEM  analysis  in  the  following  section  to  obtain  the  residual  stress  profile  in  the  20  layer  part,  the  results 
of  which  are  given  in  Table  4. 

Finite  Element  Modeling 

The  results  of  the  cutting  and  strain  gage  monitoring  summarized  in  the  previous  section  will  now  be 
combined  with  a  finite  element  modeling  (FEM)  analysis  to  obtain  the  residual  stress  profile  in  the  20  layer 
composite  part.  The  finite  element  analysis  was  completed  using  IDEAS  Master  Series  4.0  (SDRC,  1990),  a 
commercially  available  finite  element  modeling  package,  running  on  a  Hewelett  Packard  Apollo  Series  700 
workstation.  The  purpose  of  this  analysis  is  to  be  able  to  use  the  strain  gage  output  from  the  step-by-step  cutting 
process  already  completed  to  obtain  a  residual  stress  value  at  each  of  the  ten  cuts.  This  is  accomplished  by 
building  a  finite  element  model  that  is  identical  to  the  20  layer  composite  part,  including  accurate  dimensions 
and  material  properties.  Then,  unit  loads  are  applied  in  the  finite  element  model  at  the  cut  locations,  and  the 
resulting  strains  at  the  real  strain  gage  location  are  recorded  (refer  to  Figure  4).  The  strain  values  from  IDEAS  at 
the  real  strain  gage  location,  due  to  the  FEM  unit  loads,  can  then  be  correlated  with  the  known  strain  gage  output 
to  find  the  residual  stresses  in  the  real  part.  The  first  two  steps  of  this  modeling  and  analysis  procedure  are  given 
later  to  help  illustrate  the  process. 

IDEAS  20  Laver  Part  Model 

The  finite  element  model  created  in  IDEAS  is  meant  to  replicate  the  actual  part  as  closely  as  possible. 
This  model  will  be  used  in  subsequent  sections  to  obtain  the  residual  stress  profile  in  the  real  part.  The  IDEAS 
finite  element  model  for  the  20  layer  composite  part  is  as  follows: 

•  Analysis  Units  =  cm  /  N 

•  Four  Noded  Quadrilateral  Shell  Finite  Elements  (Plane  Stress) 

•  IM7/977-2  Orthotropic  Material  Properties  (given  in  the  following  section) 

•  Part  Inner  Radius  =  6.585  cm 

•  Part  Thickness  =  0.3175  cm 
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•  Part  Length  =  5.08  cm 

•  Finite  Element  Nominal  Size  =  0.03175  cm  x  0.03175  cm  square 

•  1 0  finite  elements  through  the  part  thickness 

•  Element  material  x-direction  is  oriented  along  the  arc  (circumferential) 

The  use  of  plane  stress  elements  instead  of  plane  strain  elements  will  lead  to  an  error  in  the  FEM  results 
of  less  than  1%  due  to  the  orthotropic  nature  of  the  IM7/977-2  material  (i.e.  l/O-v^-v^)  =  1.00532  =  1.0) 
(Boresi,  1993).  This  small  of  an  error  is  not  significant  for  this  analysis.  In  addition,  the  error  will  be  consistent 
throughout  the  FEM  analysis.  Thus,  because  of  the  ratio  method  used  in  the  FEM  analysis  (Equation  2),  this 
will  not  lead  to  an  error  in  the  FEM/  step-by-step  cutting  methodology  results. 

Finite  Element  Analysis  IM7/977-2  Material  Properties 

The  relevant  IM7/977-2  orthotropic  material  properties  used  in  the  IDEAS  finite  element  analysis  are 
given  below  (ICI  Fiberite,  1995).  The  x-direction  allowable  stresses  are  given  to  illustrate  that  the  residual 
stresses  in  the  part  (Table  4)  are  well  below  the  failure  values.  As  such,  the  linear  elastic  assumption  is  valid  for 
this  analysis.  The  material  coordinate  axes  (x,  y,  z)  are  set  up  in  the  FEM  so  that  the  x-direction  is  along  the 
fiber  (the  circumferential  or  hoop  direction),  the  y-direction  is  perpendicular  to  the  fiber  (the  axial  direction), 
and  the  z-direction  is  also  perpendicular  to  the  fiber  (the  radial  direction). 

•  Modulus  in  fiber  direction,  Ex  =  1.72369xl05  MPa 

•  Modulus  perpendicular  to  fiber,  Ey  =  Ez  -  1.01353xl04  MPa 

•  Poisson’s  ratio,  vxy  =  0.30 

•  Poisson’s  ratio,  vy2  =  vxz=  1.764xl0'2 

•  Shear  modulus,  Gxy  =  Gyz  =  Gxz  =  6.27423 xlO3  MPa 

•  Allowable  stress  in  tension,  x-direction  =  2.81996xl03  MPa 

•  Allowable  stress  in  compression,  x-direction  =  1.61337xl03  MPa 

A  picture  showing  the  finite  element  model  of  the  whole  part  is  given  in  Figure  2  on  the  following 
page.  The  arrows  at  the  bottom  represent  the  fixed  condition  of  the  part  along  the  bottom.  The  grid  is  the  finite 
element  mesh  of  the  part.  The  line  at  the  center  of  the  part  is  used  to  mark  the  location  of  the  centerline  of  the 
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cut.  This  line  (i.e.  the  cut)  is  approximately  2.54  cm  from  both  the  fixed  and  free  end  as  in  the  actual  20  layer 
test  part. 


Figure  2  -  IDEAS  20  Layer  Part  Finite  Element  Model 

An  enlarged  view  of  Figure  2  is  given  in  Figure  3  on  the  following  page.  This  shows  the  finite  element 
mesh  near  the  cut  in  better  detail.  In  addition,  in  Figure  3  the  first  row  of  elements  is  removed  and  unit  loads 
(arrows)  are  applied  to  model  the  first  cut. 
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Figure  3  -  IDEAS  20  Layer  Part  Finite  Element  Model  / 

Close  Up  of  Cut  Region  with  First  Cut  Modeled 

In  Figure  3  above,  one  column  by  two  rows  of  elements  (i.e.  two  elements)  are  removed  on  the  inner 
face  of  the  model  to  reflect  the  first  cut  of  the  step-by-step  cutting  procedure.  The  removal  of  one  element 
horizontally  reflects  material  removed  in  the  first  cut  (element  width  =  0.03175  cm  =  l/10th  of  the  part 
thickness).  In  addition,  unit  loads  are  applied  to  the  cut  faces  to  model  the  residual  stresses  in  the  real  part. 

20  Laver  Part  IDEAS  FEM  Analysis 

The  step  by  step  procedure  used  to  correlate  the  strain  readings  from  the  cutting  procedure  and  the 
finite  element  analysis  is  given  below.  This  procedure  is  best  explained  by  example.  Therefore,  the  first  two 
steps  of  the  actual  analysis  for  the  20  layer  part  are  given  below,  followed  by  the  general  procedure  and 
equations  that  govern  the  rest  of  the  analysis. 

After  accurately  modeling  the  20  layer  part  in  IDEAS,  the  first  step  is  to  model  the  first  cut  by 
removing  elements  from  the  mesh  and  applying  unit  loads  to  the  faces  of  the  cut.  This  is  shown  in  Figure  3 
above.  This  model  is  then  run  through  the  IDEAS  post  processing  software  to  obtain  the  strains  in  the  material- 
x  (circumferential)  direction.  The  results  of  this  analysis  are  shown  in  Figure  4  on  the  following  page. 
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Figure  4  -  IDEAS  Model  of  20  Layer  Part  /  First  Cut  with  Unit  Loading 


The  FEM  strain  values  due  to  the  unit  load  applied  at  the  cut  are  given  at  the  location  of  the  real  strain 
gage  on  the  outside  of  the  part  (i.e.  opposite  the  cut)  in  Figure  4  above.  In  addition  to  this,  the  real  strain  gage 
output  obtained  from  the  C&C  cutting  experiment  is  given  in  Table  1.  We  can  use  these  two  strain  values  (FEM 
and  real)  to  obtain  the  actual  forces  (residual  stresses)  that  were  relieved  in  the  real  part  due  to  the  first  cut. 
Because  of  the  similarity  of  the  FEM  and  the  real  part,  the  unit  loads  and  the  residual  stresses  are  related  in  the 
same  proportion  as  the  FEM  strain  and  the  real  strain  gage  strain.  This  leads  to  the  following  equation: 

x,  =eSg'  /e,"  (Equation  2) 

Where: 

x,  =  residual  stress  in  the  part  at  the  1st  cut  (the  desired  quantity) 

Esc1  =  total  strain  from  the  real  strain  gage  recorded  after  the  1st  cut 
e,"  =  IDEAS  strain  due  to  the  1st  load  (unit  load)  applied  at  the  1st  cut 
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The  value  for  -  -\66A2x\0~6  is  given  in  Table  1.  The  value  for  e/1  =  -4.520xl0'9  is  obtained  by 
taking  a  simple  mathematical  average  of  the  strain  values  in  Figure  4.  The  fact  that  both  of  these  values  are 
negative  indicates  that  they  are  compressive  strains.  Positive  values  would  indicate  tensile  strains.  The  IDEAS 
strain  values  from  Figure  4  are  summarized  in  Table  2  below  to  illustrate  how  we  obtain  e/1. 


Table  2  -  IDEAS  FEM  Analysis  Data  /  Cut  1  (20  Layer  Part) 


First  Cut 

Load  1 

20lc1io.mf1 

Nodes  Right 

Mat'l-x  Strain 

0 

-5.20E-09 

1 

-5.10E-09 

2 

-4.85E-09 

3 

-4.47E-09 

4 

-4.00E-09 

5 

-3.50E-09 

Avg.  = 

-4.520E-09 

-1.661E-04 

=  Real  Gage 

Avg.  Strain 

xi  = 

253.4044  MPa 

The  X!  =  253.4044  MPa  value  given  in  Table  2  above  was  computed  by  plugging  the  appropriate  values 
into  Equation  2.  The  Xj  value  represents  the  uniform  residual  stress  that  was  relieved  in  the  part  on  the  faces  of 
the  first  cut,  due  to  the  first  cut.  Because  the  unit  load  in  IDEAS  is  acting  toward  the  cut  faces  (refer  to  Figure 
4),  which  represents  the  relieved  residual  stress,  the  original  residual  stress  in  the  part  (Xj)  is  tensile.  Thus,  a 
positive  Xj  value  indicates  a  tensile  residual  stress  in  the  part,  and  conversely,  a  negative  X;  value  indicates  a 
compressive  residual  stress. 

We  now  have  the  average  residual  stress  in  the  part  for  the  first  l/lO111  of  the  part  thickness  (x^.  It  is 
important  to  note  that  this  analysis  procedure  gives  the  average  (i.e.  uniform)  stress  value  over  the  cut  interval. 
We  assume  this  constant  value  when  we  apply  a  constant  unit  load  to  the  cuts  in  the  FEM  model.  The  actual 
profile  is  probably  not  constant.  However,  the  true  profile  should  be  able  to  be  modeled  accurately  if  the  cuts 
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Figure  6  on  the  following  page. 
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The  value  for  esg2  =  -372.82X10-6  is  given  in  Table  1.  The  values  for  Ej21  =  -1.094xl0*8  and  Ej22  =  - 
5.473x1 0*9  are  obtained  in  the  same  manner  as  in  step  1  by  taking  a  simple  mathematical  average  of  the  strain 
values  in  Figure  5  and  Figure  6  respectively.  The  IDEAS  strain  values  and  result  for  x2  are  summarized  in  Table 
3  below. 


Table  3  -  IDEAS  FEM  Analysis  Data  /  Cut  2  (20  Layer  Part) 


Second  Cut 

Load  1 

Load  2 

201c1io.mf1 

20lc1io.mf1 

Nodes  Right 

Mat'i-x  Strain 

Mat’l-x  Strain 

0 

-1.24E-08 

-6.15E-09 

1 

-1.22E-08 

-6.06E-09 

2 

-1.16E-08 

-5.80E-09 

3 

-1.08E-08 

-5.43E-09 

4 

-9.87E-09 

-4.96E-09  1 

5 

-8.79E-09 

-4.44E-09 

i 

Avg.  = 

-1.094E-08 

-5.473E-09 

i 

-3.728E-04 

=  Real  Gage  Average  Strain 

x2  = 

-37.0092  MPa 

The  above  analysis  illustrates  the  first  two  steps  of  the  methodology  that  was  developed  to  find  the 
residual  stresses  in  the  20  layer  test  part.  The  following  is  the  general  equation  that  can  be  used  in  the  analysis 
procedure: 

X,  =  (eSG‘  -  s/'-x,  -  e“-x2 -  e,i3-Xj  - ...  -  / e,“  (Equation  4) 

> 

Where: 

Xj  =  residual  stress  in  the  part  at  the  i'th  cut  (the  desired  quantity) 

8sg  =  total  strain  from  the  real  strain  gage  recorded  after  the  i'th  cut 

e^  =  IDEAS  strain  due  to  the  j’th  unit  load  applied  after  the  i'th  cut 

In  Equation  4  above,  the  load  (j)  is  always  numbered  starting  from  1  at  the  cut  face  of  the  part  to  (i)  at 
cut  number  (i).  The  analysis  proceeds  step  by  step  through  the  piece,  using  the  previous  (x{)  values  to  obtain  the 
next.  For  this  reason,  the  experimental  error  will  accumulate  as  one  proceeds  through  the  part.  This  points  out  a 
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balance  that  must  be  used  in  the  analysis  procedure.  More  cuts  (i.e.  steps  in  the  analysis  procedure)  will  lead  to 
a  greater  accumulated  error.  However,  fewer  cuts  will  take  away  from  the  desired  resolution  in  the  residual 
stress  profile.  These  two  conflicting  influences  were  taken  into  account  and  an  analysis  procedure  was 
developed  to  minimize  the  accumulated  experimental  error  in  the  results  while  maximizing  the  resolution  of  the 
residual  stress  profile.  This  procedure  is  given  by  the  following:  The  methodology  of  Equation  4  is  used  to 
determine  the  residual  stress  values  (x;’s)  in  the  20  layer  test  part  for  the  first  five  cuts  from  the  outside-in.  This 
process  is  then  repeated  for  the  first  five  cuts  from  the  inside-out.  When  combined,  this  will  give  the  complete 
residual  stress  profile  through  the  thickness  of  the  20  layer  part  while  minimizing  the  error  in  the  results  and 
maximizing  the  resolution  of  the  residual  stress  profile. 

This  combined  methodology  was  also  used,  instead  of  analyzing  the  part  through  the  full  thickness 
from  inside  to  outside,  because  of  the  reduction  in  finite  element  modeling  required.  For  each  cut  (n)  added  to 
the  process,  (n)  additional  FEM  loading  conditions  must  be  analyzed.  The  results  of  the  IDEAS  residual  stress 
profile  analysis  for  the  20  layer  part  are  presented  in  the  following  section. 

IDEAS  FEM  Results  (20  Laver  Part) 

The  results  of  the  IDEAS  finite  element  modeling  procedure  for  the  20  layer  part  are  given  in  Table  4 
below.  The  table  reflects  the  residual  stresses  (x,)  at  each  cut  labeled  as  1  to  10  from  the  inside  to  the  outside  of 
the  part.  As  was  mentioned  previously,  a  positive  xi  value  indicates  a  tensile  residual  stress  in  the  part,  and 
conversely,  a  negative  xs  value  indicates  a  compressive  residual  stress. 
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Table  4  -  Residual  Stress  Profile  From  IDEAS  Analysis  (20  Layer  Part) 


INSIDE 

1 

253.4044  MPa 

2 

-37.0092  MPa 

3 

-61.6092  MPa 

4 

-148.8423  MPa 

5 

-89.5525  MPa 

6 

-80.1752  MPa 

7 

-40.4260  MPa 

8 

10.2471  MPa 

9 

20.8638  MPa 

10 

136.9772  MPa 

OUTSIDE 

IDEAS  FEM  Results  Validation  (20  Laver  Part) 

The  values  for  the  residual  stresses  in  Table  4  above  need  to  be  checked  to  ensure  their  accuracy.  To 
do  this,  we  note  that  the  residual  stress  profile  of  our  20  layer  part  should  have  no  net  force  (EF=0)  or  net 
moment  (SM=0)  through  the  part  thickness.  This  is  due  to  the  fact  that  the  composite  part  has  been  cut  and  has 
sprung-in;  and  as  such  it  cannot  sustain  a  net  moment  or  force.  Thus,  we  can  check  the  accuracy  of  the  residual 
stress  profile  results  in  Table  4  by  determining  the  sum  of  forces  (EF)  and  sum  of  moments  (EM)  of  the  stress 
profile.  If  these  sums  are  close  to  zero,  we  can  assume  that  the  profile  of  Table  4  is  accurate.  The  results  of  this 
analysis,  along  with  the  relative  error  from  zero,  are  given  in  Table  5  below. 


Table  5  -  Residual  Stress  Profile  EF  and  EM  (20  Layer  Part) 


Residual 
Stress  (MPa) 

Arm  from 
Center- 
line  (cm) 

Moment  = 
StressxArm 
(KN/m) 

INSIDE 

253.4044 

0.142875 

362.0515 

-37.0092 

0.111125 

-41.1265 

-61.6092 

0.079375 

-48.9023 

-148.8423 

0.047625 

-70.8861 

-89.5525 

0.015875 

-14.2165 

-80.1752 

-0.015875 

12.7278 

-40.4260 

-0.047625 

19.2529 

10.2471 

-0.079375 

-8.1336 

20.8638 

-0.111125 

-23.1849 

OUTSIDE 

136.9772 

-0.142875 

-195.7062 

EF  = 

-36.1219 

MPa 

EM  = 

-8.1239 

KN/m 

EF  /  X|F|  = 

-4.11 

% 

EM  /  E|M|  = 

-1.02 

% 
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As  can  be  seen  in  the  last  row  of  Table  5,  the  error  in  the  EF  and  EM  calculations  is  relatively  small.  In 
addition,  the  more  critical  measure  for  our  analysis,  the  EM  error,  is  only  1%.  This  suggests  that  the  step  by  step 
cutting/FEM  analysis  for  the  20  layer  part  has  yielded  accurate  results  for  the  residual  stress  profile  (Table  4). 

Conclusion 

This  paper  has  outlined  a  step  by  step  cutting  /  finite  element  analysis  procedure  that  can  be  used  to 
determine  the  residual  stress  profile  in  a  thin  composite  part.  This  procedure  was  originally  developed  to  assist 
in  determining  the  spring-in  in  thin  composite  parts.  However,  this  procedure  can  easily  be  expanded  to  include 
thick  composite  parts  or  isotropic  materials. 
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A  STUDY  OF  THE  EFFECTS  OF  NOVAE  ON 
THE  INFRARED  CELESTIAL  BACKGROUND 


David  A.  Joiner 
Research  Assistant 

Department  of  Physics,  Applied  Physics,  and  Astronomy 


Abstract 

The  low  densities  and  short  timescales  involved  with  dust  grain  formation  in  nova  outflows  imply  that  the 
grain  size  distribution  of  dust  which  forms  does  not  conform  with  thermodynamic  equilibrium,  and  grain  growth 
in  these  environments  should  be  modeled  in  terms  of  kinetic  equations  (Johnson  et.  al.  1993).  One  recent  use 
of  kinetic  equations  to  model  dust  growth  in  stellar  outflows  is  presented  in  Egan  and  Leung  (1995.)  We 
present  corrections  to  the  equations  for  kinetic  growth  presented  in  Egan  and  Leung  (1995),  and  compare  the 
method  of  contracting  a  large  set  of  kinetic  equations  by  taking  moments  of  the  kinetic  equation  to  a  method  of 
binning  equations. 

To  determine  the  effect  of  novae  on  the  infrared  background,  and  particularly  data  from  the  MSX  mission, 
we  first  calculate  the  expected  rate  of  nova  observations  present  in  the  MSX  catalog.  We  find  that  the  expected 
number  of  classical  novae  for  any  given  map  in  band  A  (8.28  pm)  is  of  the  order  of  10,  and  on  the  order  of  1  in 
all  other  bands.  We  compare  1  known  observation  in  the  visual  which  brightened  shortly  before  MSX 
measurements  of  it’s  location.  We  report  a  band  A  measurement  of  .322  Jy  for  Nova  V4361  Sgr  (Sgr  1996), 
and  discuss  the  implications  of  this  measurement  relative  to  visual  measurements  published  on  this  nova. 


4-2 


Introduction 

Long  ago  when  curious  minds  looked  at  the  sky  and  saw  a  light  appear  where  before  there  had  been  only 
emptiness,  they  thought  they  were  witnessing  the  birth  of  a  new  star.  They  even  gave  it  the  name  ’’nova  Stella”, 
or  ’’new  star”  (Bode  and  Evans,  1987).  The  phenomenon  of  the  nova  is  now  known  to  be  part  of  a  broad  class 
of  stars  known  as  cataclysmic  variables,  in  which  large  changes  in  brightness  occur  in  very  brief  time  periods. 
In  the  classical  nova,  matter  accretes  from  a  binary  companion  onto  the  surface  of  a  white  dwarf  star.  The 
heating  of  this  material  by  the  white  dwarf  erupts  in  a  thermonuclear  runaway,  launching  a  shell  of  hydrogen, 
carbon,  and  other  elements  at  ~  1000km  s-',  and  the  star  is  observed  to  brighten  by  many  magnitudes  in  the 
infrared,  visible,  and  ultraviolet  (>  9  magnitudes  in  the  visible  spectrum.)  Peak  luminosities  in  the  visible  can 
range  from  104  -  106L©  for  classical  novae,  with  anywhere  from  a  few  thousandths  to  all  of  the  total 
luminosity  reradiated  in  the  infrared  (Stanfield  1988). 

Over  the  course  of  the  nova’s  evolution,  great  changes  can  be  seen  in  the  infrared,  visible,  ultraviolet,  and 
x-ray  regions  of  the  spectrum.  The  only  major  difference  between  visible  light  curves  is  the  possible  occurrence 
of  a  deep  dip  in  the  visible  light  curve  during  the  decline  period,  associated  with  an  increase  in  the  infrared. 
Both  changes  are  on  the  order  of  a  few  magnitudes.  The  cunent  accepted  explanation  for  this  dip  is  that  it  is 
caused  by  the  formation  of  dust  in  the  nova  outflow. 

The  nova  event  generally  lasts  for  ~  100  days,  and  the  infrared  excess  associated  with  dust  formation 
generally  takes  place  over  a  few  days  or  weeks.  Approximately  100  novae  occur  in  the  galaxy  every  year,  of 
which  about  20%  have  an  apparent  visible  magnitude  above  6.0.  While  the  infrared  signature  during  the  initial 
expansion  is  due  only  to  the  expanding  shell  of  gas,  the  infrared  signature  of  the  dust  formation  event  is 
dominated  by  the  absorption  and  radiation  of  carbonaceous  or  silicate  material  in  the  dust  shell.  In  order  to 
understand  the  spectral  signature  of  the  dust  formation  period  in  classical  novae,  the  growth  of  these  grains 
must  then  be  modeled,  and  the  radiative  transfer  of  light  through  this  dusty  shell  determined  from  the  dust 
characteristics. 

Solving  the  kinetic  equations  which  describe  the  growth  of  dust  grains  has  primarily  one  obstacle.  A  1  pm 
sized  grain  requires  1012  equations.  Researchers  generally  reduce  this  set  of  equations  by  one  of  two  methods, 
taking  moments  of  the  equations,  or  grouping  equations  into  bins. 

As  part  of  this  study  I  have  also  estimated  the  number  of  novae  present  in  the  MSX  satellite  data  is 
estimated.  Visual  observations  of  the  slow  nova  V4361  Sgr  are  interpreted  with  respect  to  an  infrared 
observation  (8.28pm)from  the  MSX  point  source  catalog.  The  density  of  the  shell  on  July  19,  1996,  is 
calculated  to  be  ~  10,0cm-3. 

Modeling  Grain  formation 
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Modeling  the  kinetics  involved  in  grain  formation  require  the  solution  of  a  large  set  of  nonlinear  ordinary 
differential  equations.  The  general  approach  is  to  in  some  way  contract  the  large  numbers  of  equations  required 
to  a  number  which  is  within  one’s  computational  means,  and  integrate  using  a  method  suitable  for  stiff  sets  of 
equations,  such  as  the  backwards  Euler  method  (Press  et.  al.  1992)  or  the  Gear  method  (Kahaner  et.  al.  1989). 
Improvement  of  computer  code  from  Egan  and  Leung. 

Egan  and  Leung  (1995)  cite  the  following  equations  for  the  kinetic  growth  of  grains,  determined  via  a 
method  of  taking  moments  of  the  kinetic  equations,  coupled  to  a  truncated  set  of  equations  up  to  some  cutoff: 


f  =  -  Pu-jfi  +  Pi+j.ifi+j  +  Pi-ufi-j]  '  > 


M 


f  =  HUji +  Pjj-ifj]  i  *  M 

+  i-m  ~  P i  ^ i+m,ifi i+m  ~~  P i,i+mf  f  ] 

dKj  | 
dt 


where  Pu+j  is  the  reaction  rate  for  the  kinetic  equation  Cj  +  C,  -►  Cy+/,  fj  is  the  density  of  the  C}  gas 
component,  and  only  growth  by  accretion  of  clusters  up  to  size  M  is  considered.  I  is  a  factor  which  takes  into 
account  the  difference  in  the  rate  for  the  collision  between  two  clusters  of  the  same  size,  and  as  it  is  written  in 
the  equations  above  is  given  by  l  =  (1  +  Suj)-  H‘  is  the  interaction  of  grains  above  the  cutoff  grain  size  Nc 
with  the  ith  grain  size.  In  the  formalism  presented  by  Egan  and  Leung  for  determining  -^L,  each  i  term  is 
linearly  independent,  and  the  determination  of  -!,•  from  is  trivial. 

This  expression  for  the  accretable  species  neglects  reactions  of  the  type  C,-  +  Cm  ^  Cj+m,  and  also  has 
double  counting  for  reactions  with  grains  larger  than  Nc .  If  the  total  shell  mass  is  given  by  M  -  V^^={  mtfi, 
and  assuming  that  the  volume  V  is  constant,  and  the  mass  of  the  ith  grain  size  is  mf-  =  imc ,  then  conservation  of 
mass  requires  that 


oo  r  Nc~l  co  — 

i-Lv,  =  i 

i=l  L  *=1  i—Nc  J 


Nc-\ 


X) ift  +^3  =0 


1=1 


Taking  these  factors  into  consideration,  the  improved  set  of  kinetic  equations  is  given  by 


Ml 

dt 


|  i,i+jfi  P i  +  Pi+j\if i+j  PH,if  H  ] 


/  >  M 
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-  HjUJi\rpHdfH  + 

+  j  fm  [P f-m  “  P i,i-mf  r  +  Pi+m,if i+m  P i,i+mf f  ] 
+  i+m, if i f+m  —  ^  U+mf i i  ] 

+t*  Ej,>’ -  +  ^/v+t^] 

dKy  , 

”  3? 


i  <  M 


k>  - 

These  alterations  were  made  in  Egan’s  time  dependant  grain  growth  code,  and  are  currently  being  tested. 

Comparison  of  bin  and  moment  method 

As  was  described  in  the  previous  section,  the  kinetic  growth  of  dust  grains  by  accretion  of  small  clusters 
has  been  done  primarily  by  two  methods.  This  was  originally  attempted  by  Yamamoto  and  Nishida  (1977),  but 
they  did  not  use  any  method  to  contract  their  system  of  equations,  and  their  grain  size  was  limited  to  the  number 
of  equations  they  could  solve  practically,  which  with  the  computing  resources  at  that  time  was  ~  100  equations. 
Gail  and  Sedlmayr  (1988)  introduced  a  closed  set  of  moments  of  the  kinetic  equations.  Egan  and  Leung  (1995) 
applied  this  method  to  the  growth  of  dust  grains  in  the  outflows  of  Asymptotic  Giant  Branch  stars. 

A  second  approach  was  taken  by  Johnson,  Friedlander,  and  Katz  (1993),  in  which  the  equations  are 
grouped  into  bins  when  the  grain  size  becomes  large. 

For  the  purposes  of  this  comparison,  growth  by  accretion  of  monomers  only  will  be  considered  for 
simplicity. 

Consider  the  following  model 

Ei  =  (1  -0.5$f-i,i)4jr(/-i-i  +ao)2\~ — % — TtI 

|_  {znmckl  grain)  J 

Ci  =  (l-0.55u)*(n  +  «o)2(-^)l/2 

where  C,-  is  Pu+\  and  E,  is  Py-ifi  in  our  previous  notation,  -y-  is  a  normalization  factor  allowing  us  to  write 
the  equation  in  terms  of  total  numbers  instead  of  densities.  This  gives  us  the  following  master  equations  to  be 
solved 

■*f-=  N ,  [-C/JV,  +  Cj-\Ni-\ ]  i  >  1 

+[-EiNj  +  Ej+iNj+i] 

dN  1  _  _  y-'00  -•  dN,  •  „  1 

dt  -  dt 


Moments  can  be  taken  of  this  set  of  equations  to  reduce  the  number  of  equations  needed  to  be  solved.  For 
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the  above  definitions,  if  i  >  Nc  then  Ni  »  Ni+ 1, 


M  [ — C/7V +  Cj-\Ni~\  ]  +  +  EmNm ]  j  >  1 

+[-£'l7V(-  +  Ei+\Ni+i  ] 


dN  ] 
dt 


A/=2 


I  AO-  -*3 


-fL= 


</*,- 

</r 


1/2 


/Vo 


Vi  +  4^ao 


0  (2nmckT grain)1?2 


i  >  0 


Alternatively,  if  the  master  equations  are  expressed  in  binned  form  with  i  grain  sizes  g,-,  then  the  rate  of 
growth  from  bin  i  to  bin  i  +  1  is  given  by  the  rate  of  accretion  of  (g,+i  -  g/)  monomers  assuming  that  the  rates 
are  essentially  unchanged  until  growth  to  the  next  bin  has  occurred.  The  binned  equations  can  be  written  as 


dN, 

dt 


w,[ 

+[" 


Ci  Ni  cMjyM 

8i+\~8i  gi-gi-l 

EjNj  EmNm  “I 

8i~8i- 1  Si+l^gi  J 


] 


i  >  1 


dN  i  _  dN, 

dt  “  Zjj=2°'  dt 


i  =  1 


The  above  problem  is  solved  using  each  method.  Since  the  moment  method  does  not  explicitly  give  us  the 
grain  size  distribution,  moments  are  taken  of  the  grain  size  distribution  from  the  bin  method,  and  these  are 
compared. 

The  parameters  for  the  comparison  models  are  given  in  table  1 . 


ngas(r  =  rt) 

2.893  x  106cm~3 

U 

1O4L0 

Tt 

2500  K 

T gas  =  T grain 

r*  (/?*//?) 1/2 

v„u, 

106 

Table  1. 


Values  typical  of  an  AGB  star  are  used  for  comparison  with  the  work  by  Egan  and  Leung.  All  integrations 
were  performed  using  Kahaner  and  Sutherland’s  DDRTV2  subroutine,  in  Fortran77,  on  a  Dec  Alpha  with  a  225 
MHz  processor  (Kahaner  et.  al.  1989).  The  time  and  memory  requirements  of  each  method  depend  primarily 
on  the  total  number  of  equations  required  in  the  solution. 

Results  for  each  method  agreed  qualitatively,  and  Figure  1  shows  grain  size  as  a  function  of  radius  from  the 


4-6 


central  source. 

The  major  approximation  that  is  made  in  the  derivation  of  the  moment  equations  is  the  assumption  that  for 
large  i,  Nj  ~  Nj-i.  As  Nc  is  increased,  this  approximation  should  improve,  and  the  integration  result  should 
approach  the  exact  solution. 

In  the  bin  method,  the  approximation  is  that  a  discrete  set  of  grain  sizes  is  treated  as  a  set  of  bins  above 
some  cutoff  Nb.  The  grid  spacing  after  that  is  arbitrary,  but  we  have  chosen  a  logarithmic  bin  spacing.  The 
solution  should  approach  the  exact  solution  as  the  bin  widths  decrease.  There  are  then  two  parameters  which 
can  be  changed  to  improve  the  accuracy  of  the  results:  the  value  of  the  cutoff  at  which  binning  begins,  and  the 
logarithmic  spacing  factor  x- 
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Figure  1 .  Grain  size  as  a  function  of  radius  for  the  bin  method  with  Nb  =  90,  75 
logarithmically  spaced  equations,  and  a  maximum  possible  grain  size  of  108. 

Figure  2  shows  the  maximum  grain  size  from  each  model,  plotted  against  the  total  number  of  equations 
required  in  the  solution  for  4  cases.  Case  1  is  the  moment  method.  Case  2  is  the  bin  method,  with  the  cutoff 
value  N/,  held  constant.  Cases  3  and  4  show  how  the  solution  changes  with  N*,  if  the  logarithmic  spacing  factor 
X  is  held  constant.  The  total  number  of  equations  required  in  the  moment  method  depends  on  the  number  of 
extra  kinetic  equations  above  Nc  used,  and  on  the  number  of  moments  solved  for.  For  these  models, 
N moment  =  Nc  +  7.  In  the  bin  method,  the  number  of  equations  beyond  Nb  is  determined  by  and  is  set  so  that 
the  size  of  grains  in  the  last  bin  is  108  carbon  atoms. 
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All  methods  show  an  increase  in  grain  size  with  an  increase  in  the  accuracy  of  the  equations.  The  bin 
method  seems  to  be  leveling  off  at  a  smaller  number  of  equations.  It  is  not  yet  clear  exactly  if  the  moment 
method  will  converge  to  the  same  solution,  and  this  will  require  further  calculations  to  determine.  If  we  assume 
that  all  models  are  approaching  the  same  value,  which  due  to  the  monotonic  increase  in  each  curve  must  be 
greater  than  the  values  plotted,  then  it  appears  that  the  bin  method  is  capable  of  achieving  greater  accuracy  for 
the  same  number  of  equations  used. 


Figure  2.  Average  grain  size  plotted  against  number  of  equations  used 
for  moment  and  bin  methods.  For  the  bin  method,  either  the  logarithmic 
spacing  for  the  binned  equations  or  the  grain  size  at  which  binning  begins 

is  held  constant 

Novae  in  Celestial  Backgrounds 

The  continuous  novae  spectrum  can  generally  be  described  in  three  phases.  At  outburst,  an  optically  thick 
fireball  expands  at  constant  velocity  and  constant  bolometric  luminosity,  with  a  decreasing  effective 
temperature.  As  the  shell  becomes  optically  thin,  the  continuous  spectra  consists  of  an  optically  thin  blackbody 
in  the  visible,  free-free  excess  at  optically  thin  wavelengths,  and  an  optically  thick  blackbody  at  long 
wavelengths,  where  the  optical  depth  is  dominated  by  free-free  absorption.  The  effective  temperature  of  the 
nova  rises,  while  the  bolometric  luminosity  remains  constant,  and  the  visible  flux  decreases.  For  many  novae, 
particularly  those  with  higher  mass  loss  and  lower  luminosity,  a  dust  formation  event  takes  place  which 
produces  an  optically  thick  shell  in  the  visible,  radiating  as  a  ~  1 500  K  blackbody  in  the  infrared. 

Parameters  for  a  typical  nova  that  will  be  used  in  determining  the  rate  of  occurrence  of  novae  in  the  MSX 
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data  are  listed  in  table  2. 

For  the  temperature  of  the  novae,  T  is  determined  by  the  radius-luminosity-temperature  relation  when  the 
shell  is  optically  thick,  and  is  picked  so  that  the  decline  of  the  visible  luminosity  matches  a  decline  speed  of 
r2  -  50  days  in  the  optically  thin  phase. 


Outflow  velocity 

Vou, 

5.0  x  107cm  s'1 

Sound  speed 

V sound 

1.5  x  106cm  s-! 

Shell  radius 

R 

V OUtt 

Shell  thickness 

l 

V sound? 

Luminosity 

L 

10  5Le 

Effective  Temperature 

T 

z  <  1  To{ULQ)w{RIRQ)-'a 
z  >  1  determined  by  f2  ~  50  days 

Grain  Temperature 

h 

T{RJR)m 

Mass  Loss 

M 

10  -5M© 

Carbon  Mass 

Mc 

5  x  10  ~2M 

Grain  size 

a 

0.3  pm 

Table  2. 


The  flux  is  modeled  as 

Fx  =  (l-e-r)Bx*A  +  eF/  *V*e~r 

where  r  is  the  optical  depth  at  the  given  wavelength,  A  is  the  area  of  the  shell,  V  is  the  volume  of  the  shell,  and 
e£F  is  the  free-free  emissivity  (Rybicki  and  Lightman  1979), 

CFF  dv  rFF  -  c  (  25ne6  r  2 n  \  1/27-1/272-  „  r-hvfkT^  \ 

^  "l^(l£r(31bi«)  r  Zn‘n'e  **> 

Approximately  half  of  all  novae  form  dust  shells,  so  for  the  purposes  of  these  calculations  a  typical 
classical  novae  has  a  50%  chance  of  forming  an  optically  thick  dust  shell  at  Tg  =  1500 K.  It  is  also  assumed 
that  the  hydrogen  is  completely  ionized. 

Figures  3  and  4  show  the  average  flux  in  each  band  for  MSX  bands  A,  Bl,  B2,  C,  D,  and  E.  The  values 
plotted  are  D2FV,  where  D  is  the  distance  to  the  source  in  centimeters,  and  Fv  is  the  flux  in  Janskys.  Figure  3 
applies  to  a  nova  with  an  optically  thick  dust  shell  forms,  and  figure  4  applies  to  a  nova  which  does  not  form 
dust. 

The  average  rate  of  novae  in  our  galaxy  is  lOOy*"1  (Bode  and  Evans  1987).  If  the  majority  of  these  novae 
occur  in  the  plane,  and  the  solid  angle  subtended  by  the  plane  of  the  galaxy  is  approximately  0.2  rad  by  2k  rad, 
or  -  1 .25  sr,  this  yields  a  rate  per  year  per  steradian  in  the  galactic  plane  of  80  y-1  sr1.  Assuming  these  novae 
are  evenly  distributed  along  a  line  of  sight  out  to  20  kpc,  but  the  area  subtended  by  a  solid  angle  passing  through 
a  flat  cylinder  is  proportional  to  R,  then  the  probability  function  of  a  novae  occurring  spatially  will  be 
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p(r )  =  (5  x  10-3kpc-2)/?,  where  integration  of  p(r)  out  to  20  kpc  has  been  normalized  to  1.  The  rate  of 
occurrence  of  novae  at  a  given  radius  looking  out  along  a  line  of  sight  is  then  r^1  =  0.4  y-1  sr_1kpc ~2RdR.  The 
timescale  for  these  events  will  be  days,  which  is  much  longer  than  the  duration  of  the  measurements,  so  the 
number  actually  seen  will  be  the  rate  of  occurrence  multiplied  by  the  duration  of  the  event. 


t(s) 


Figure  3.  Band  A  flux,  normalized  by  the  distance  to  the 
source,  for  a  typical  dust  producing  nova. 
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Figure  4.  Band  A  flux,  normalized  by  the  distance  to  the 
source,  for  a  typical  non-dust  producing  nova. 

The  duration  (Af)  that  a  novae  will  be  intense  enough  to  be  seen  by  the  instrument  will  change  as  the  novae 
occurs  farther  away  from  us,  so  the  average  number  seen  per  steradian  in  an  observation  will  be 

iV  =  0.4  J  At(R)RdR 

The  results  for  each  band  are  given  in  table  3. 


Band 

A 

B 1 

B2 

C 

D 

E 

A(sr~l) 

13.1 

.47 

.78 

1.40 

1.39 

.46 

Table  3. 


At  best,  this  should  be  considered  an  order  of  magnitude  calculation,  however,  due  to  the  assumptions 
made  throughout  the  calculation.  To  an  order  of  magnitude,  the  rate  of  nova  per  observation  per  steradian  in 
the  MSX  data  should  be  ~  10  in  band  A,  and  ~  1  in  each  of  the  other  bands. 

Nova  V4361  Sgr 


Figure  5.  Light  curve  in  V  band  for  V4361  Sgr.  All  values  are  from  LAU 
circulars,  and  authors  for  each  data  are  listed  in  table  4. 

In  the  interest  of  finding  novae  in  the  MSX  data  set,  any  nova  observed  at  visible  wavelengths  during  the 

time  of  the  MSX  mission  were  considered.  Two  novae  occurred  during  the  MSX  mission  that  were  observed  in 

the  visible,  and  occurred  in  the  galactic  plane,  but  only  one  of  these,  V4361  Sgr,  occurred  shortly  before  MSX 
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observed  that  part  of  the  sky.  V4361  was  first  observed  after  outburst  on  July  11,  1996,  at  b  =  -2.28,  /  =  13.7. 
MSX  measured  this  part  of  the  sky  once  on  July  19th,  and  later  3  times  in  early  September.  There  was  a  point 
source  detected  in  band  A  (8.28  pm),  and  no  other  bands,  in  the  July  19th  measurement,  and  no  detections  at  all 
in  the  September  data. 

Table  4  gives  data  published  in  IAU  circulars  6443,  6447,  6450,  6453,  and  6454  by  various  authors.  The 
visual  light  curve  is  plotted  in  figure  5. 


Source 

Day  (Aug  96) 

U 

B 

m 

□n 

H 

m 

H 

K 

Nakano  et.  al. 

EH 

Hi 

Gilmore  and  Kilmartin 

-17.502 

11.89 

11.69 

Pel 

Monard 

5.86 

Hi 

Morrison  and  Argyle 

5.94 

11.61 

H 

Ripero  Osorio 

6.88 

11 

■ 

Morrison  and  Argyle 

6.94 

11.2 

Dillon 

7.22 

10.8 

Schmeer 

7.9 

11.3 

Morrison  and  Argyle 

7.94 

11.54 

Morrison  and  Argyle 

8.93 

11.65 

Dillon 

9.15 

11.2 

Sostero  et.  al. 

9.924 

12.78 

11.78 

10.92 

10.14 

Morrison  and  Argyle 

9.93 

11.96 

Martin 

9.965 

8.93 

8.73 

8.35 

Sostero  et.  al. 

10.877 

13.03 

12.04 

11.19 

10.39 

Morrison  and  Argyle 

10.93 

12.29 

Morrison  and  Argyle 

11.92 

12.07 

Sostero  et.  al. 

12.855 

12.99 

12.06 

11.24 

10.45 

Morrison  and  Argyle 

12.92 

12.19 

Sostero  et.  al. 

13.917 

12.97 

11.97 

11.17 

10.38 

Morrison  and  Argyle 

13.92 

Sostero  et.  al. 

14.973 

Sostero  et.  al. 

15.868 

umi 

m9i 

Sostero  et.  al. 

16.854 

12.87 

11.81 

11.02 

10.22 

Table  4. 


Dust  formation  phase 

A  classical  novae  reaches  a  maximum  in  infrared  wavelengths  at  the  time  the  nova  shell  becomes  optically 
thin  in  thermal  Bremsstrahlung  absorption.  A  second  infrared  emission  phase  can  occur  in  novae  which  form 
dust.  V4361  did  have  a  significant  dust  formation  event.  There  were  3  MSX  maps  that  covered  the  area  of  the 
sky  in  which  V4361  Sgr  occurred  after  the  dust  formation  event.  The  optical  depth  of  the  dust  formation  event, 
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computed  by  comparing  the  minimum  of  the  light  curve  during  the  event  with  an  extrapolation  of  the  light  curve 
from  before  the  dust  formation  event,  gives  us  an  optical  depth  of  rv  =  1.1  for  the  dust  shell  at  the  time  of 
maximum  extinction.  However,  the  dust  shell  had  expanded  to  the  point  of  being  optically  thin  well  before  the 
MSX  measurements,  so  it  is  not  surprising  that  no  detection  was  made.  While  a  direct  estimate  without 
knowing  the  dynamics  of  the  shell  is  difficult,  it  is  known  that  the  light  curves  of  all  nova  are  similar  on  a 
dimensionless  time  scale,  i.e.,  the  ratio  of  the  time  at  which  an  event  occurs  in  the  light  curve  to  the  outburst 
time  is  similar  for  all  novae.  If  we  consider  an  outburst  time  of  July  10th,  and  a  dust  formation  time  of  August 
8th,  and  scale  our  typical  A  band  light  curve  from  figure  3  accordingly,  we  see  that  the  flux  around  September 
10th  is  4  orders  of  magnitude  less  than  the  maximum  Bremsstrahlung  emission. 

Bremsstrahlung  emission  phase 

The  July  measurement  did  get  a  glimpse  of  emission  in  band  A,  at  A  =  8.28  pm  (AA  =  4  pm),  of  .322  Jy, 
well  above  the  sensitivity  limit  of  the  instrument  (0.07  Jy).  In  units  of  Ex,  this  is  1 .41  x  10~15  erg  cm-2  s"1  A-1 . 
Due  to  the  time  of  this  measurement  with  respect  to  the  novae  outburst,  and  the  fact  that  the  value  at  this 
wavelength  exceeds  the  blackbody  radiation  at  the  novae  temperature,  this  is  most  likely  thermal 
bremsstrahlung  emission. 

Color  color  correction 

The  nearest  spectra  to  this  date  are  the  UyByV,RJ  spectra  taken  by  Gilmore  and  Kilmarten  on  July  13.498. 
Given  the  very  slow  speed  class  {ti  ~  200 d)  of  the  novae,  and  due  to  the  fact  that  the  dust  formation  does  not 
occur  for  another  16  days,  the  spectra  at  the  time  of  the  MSX  measurement  should  be  approximately  the  same 
as  the  spectra  taken  6  days  earlier.  The  general  trend  in  novae  is  towards  a  slowly  rising  temperature,  with 
constant  bolometric  luminosity,  so  any  error  in  doing  this  would  lead  to  an  underestimation  of  the  IR  excess  due 
to  free-free  emission. 

This  data  must  first  be  corrected  for  reddening  by  foreground  material.  A  value  of  E(V-  [12])  for  the 
galactic  plane  is  taken  from  Egan  and  Price  (1996). 


£(V-[  12])  =W)  = 


0.1 


sin(lfrl)  +  0.03(2  -  cos(/)) 

Values  from  Whittet  are  used  for  the  average  interstellar  extinction  curve  (1992).  The  interstellar  extinction 
curve  follows  a  power  law  in  the  infrared,  given  by 


E(  A-V) 
E(B-V) 


1.19A”1 84  -  3.05 


holds  at  least  up  to  5  microns.  Assuming  that  the  same  power  law  holds  beyond  5  microns,  then 


E(B-V) 


iiy _ 

1 . 1 9[  1 2]_l  84  -  3.05 


329R(b,l) 
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The  position  of  V4361  Sgr  is  b=-2.28, 1=13.7,  which  leads  to 

E(B-V)  =  3.29/?(-2.28, 13.7)  =  0.466 


Band 

nu 

E(X- V) 

Ma 

fA 

?BB 

U 

11.89 

.746 

11.14 

1.49e-13 

1.35e-13 

B 

11.69 

.466 

11.22 

2.1  le-13 

1.88e-13 

V 

10.59 

0 

10.59 

2.1  le-13 

2.1  le-13 

R 

9.79 

-0.363 

10.15 

1 .5  le-13 

1.83e-13 

/ 

9.04 

-.746 

9.786 

1.01e-13 

1.27e-13 

4 

6.91  | 

-1.41 

8.32 

3.83e-16 

9.54e-17 

Table  5. 


While  exact  values  require  knowing  the  extinction  Ay,  relative  values  can  be  determined  by  using  a  value 
of  A  v  =  0.  The  blackbody  spectrum  shown  is  the  best  fit  to  the  U,B,V,  R  and  /  measurements,  using  the  method 
of  simulated  annealing,  which  matched  a  temperature  of  5330  K.  Dereddened  values  are  listed  in  table  5. 
Hydrogen  gas  density 

If  all  of  the  radiation  in  band  A  is  from  Bremsstrahlung  emission,  the  gas  density  at  the  time  of 
measurement  can  be  estimated.  The  main  spectral  line  seen  in  novae  outburst  near  band  A  is  the  NVI  line  at 
7.6  pm,  but  that  is  generally  seen  much  later  in  the  nova  outburst  (Gehrz  1988).  At  the  time  the  measurement  is 
taken,  the  ratio  of  the  infrared  emission  to  a  blackbody  extrapolated  from  a  fit  to  the  (7,B,V,/f,and  I 
measurements  yields 

Eff  E_i  _  4  m 

Ebb  Ebb 


Using  the  previous  definitions  of  the  nova  flux, 

Eff  m  (1  -  e~Tv)5y  +  SyFLg~r' 
Ebb  tvisBv 


This  function  is  solved  numerically  for  Z  =  1,7'=  5330,  and  ne  =  —  n,  to  obtain  a  relationship  between  the 

shell  thickness,  L,  and  the  gas  density,  n.  The  results  are  shown  in  figure  6. 

The  thickness  of  the  shell  needs  to  be  less  than  10 12  cm  for  the  shell  to  be  in  the  optically  thin  regime. 
Clearly  for  any  realistic  values  of  L,  the  band  A  measurement  is  taken  at  an  optically  thick  wavelength,  and 


Eff  _ 
Ebb 


-l 


For  V4361  Sgr,  ne  =  (4.0 otE)~1.  While  the  geometry  of  this  particular  novae  is  not  well  known,  given  a 
typical  shell  thickness  of  L=1013,  a  gas  density  of  4.2  x  1010  is  calculated.  Given  uncertainty  in  the  thickness 
of  the  nova  shell,  this  should  be  considered  at  best  an  order  of  magnitude,  and  n  ~  10 10  cm-3. 
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Figure  6.  Results  of  numerical  solution  for  shell  density  as 
a  function  of  shell  thickness  for  V4361  Sgr. 

Conclusions 

Grain  growth 

In  order  to  conserve  mass,  alterations  are  required  in  the  method  of  Egan  and  Leung.  We  have 
implemented,  but  are  still  in  the  process  of  testing  the  results  of  these  changes,  both  for  the  past  application  of 
this  method  to  AGB  stars,  and  to  the  current  application  of  this  method  to  the  study  of  grain  growth  in  nova 
shells. 

For  the  case  of  monomer  accretion,  while  there  is  good  qualitative  agreement  between  the  bin  and  moment 
methods,  there  can  be  as  much  as  50%  difference  in  the  results  of  calculations  for  parameters  typical  of  AGB 
stars.  We  intend  to  continue  this  testing  as  it  applies  to  novae.  Current  results  seem  to  imply  that  the  bin 
method  attains  greater  accuracy  for  a  smaller  number  of  equations.  It  should  be  noted,  however,  that  this  testing 
was  for  the  case  of  monomer  only  accretion.  It  is  not  yet  clear  how  well  suited  the  bin  method  is  to  the  growth 
of  grains  by  accretion  of  larger  clusters,  and  comparison  of  this  to  the  moment  method  would  help  to  further 
differentiate  between  these  two  methods.  Also,  while  the  bin  method  gives  explicit  information  about  the  grain 
size  distribution,  the  moment  method  requires  the  use  of  a  reconstruction  method  to  determine  the  grain  size 
distribution.  Future  work  should  compare  a  maximum  entropy  reconstruction  as  in  Vicanek  and  Ghoniem 
(1992)  with  the  explicit  distribution  from  the  bin  method. 

Novae  in  the  MSX  data 

Classical  novae  should  be  expected  to  appear  in  the  MSX  data  set,  but  not  in  great  abundance.  However, 
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the  timing  of  the  experiment  was  sufficient  to  obtain  useful  infrared  data  for  at  least  one  nova,  V4361  Sgr. 
While  the  likelihood  of  using  MSX  data  to  account  for  other  classical  novae  is  not  high,  one  might  consider  the 
effect  of  dwarf  novae  on  the  MSX  measurements,  as  their  frequency  is  higher,  and  their  eruptions  are  recurrent. 
Nova  V4361  Sgr 

We  confirm  the  sighting  of  a  known  nova,  V4361  Sgr,  in  the  MSX  data,  and  compare  the  infrared  flux  with 
other  visual  measurements  to  calculate  a  hydrogen  shell  density  of  1 0 10  cm  ®  at  July  19,  1996.  Wagner  et.  al. 
report  taking  spectra  of  hydrogen  Balmer  lines,  Fe  H,  and  Ca  H,  but  have  not  yet  published  an  analysis  of  these 
lines  to  determine  the  outflow  velocity  of  the  source  (Nakano  et.  al.  1996).  Knowledge  of  the  outflow  velocity 
will  help  to  better  constrain  these  calculations. 
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Abstract 


Blumleins  are  transmission-line  structures  that  allow  the  formation  and  propagation  of 
electromagnetic  waves.  They  have  existed  for  decades  and  their  use  has  been  widely  documented  for  a 
variety  of  applications.  The  most  common  form  of  the  single  Blumlein  is,  perhaps,  the  three-conductor 
coaxial  line.  The  widespread  use  and  the  coaxial  cable’s  simple  geometry  lend  itself  well  to  today’s  needs. 
The  tri-plate  Blumlein  configuration,  on  the  other  hand,  is  not  as  well  known.  Due  to  its  more  complex 
geometry,  certain  aspects  must  be  considered  which  are  not  at  all  relevant  in  the  coaxial  case.  This  work 
attempts  to  explore  the  intricacies  of  not  one  triplate  Blumlein,  but  a  stack  of  such  devices.  A  model  has 
been  constructed  which  not  only  accounts  for  wave  propagation  in  the  time  domain,  but  also  Blumlein 
charging  and  commutation.  Each  portion  of  the  model  has  been  compared  to  modified  (existing)  analyses 
of  similar  structures.  The  limits  of  validity  for  this  analysis  have  also  been  tested  through  experimental 
studies. 


A  COMPUTATIONAL  ANALYSIS  OF  STACKED  BLUMLEINS 
USED  IN  PULSED  POWER  DEVICES 

Johnelle  L.  Korioth 


Introduction 

Modem  applications  for  electromagnetic  waves  have  a  basis  with  the  discoveries  of  Maxwell, 
Faraday,  Ampere  and  Hertz.  There  has  been  much  success  in  most  all  frequency  ranges  with  regard  to 
fundamental  understanding  and  technological  implementation.  To  some  extent,  though,  there  still  exists  a 
gap  in  the  millimeter  and  sub-millimeter  regions.1  This  zone  is  of  great  interest  because  some  of  the  typical 
resonant  frequencies  of  many  biological  and  physical  objects  lie  in  this  domain  and,  more  simply,  because 
this  is  the  area  of  crossover  from  radio  to  optical  frequencies.  It  is  imperative  that  this  gap  be  mastered  if 
applications  like  enhanced  radar  systems,  atmospheric  chemical  analysis,  industrial  materials  processing, 
and  advanced  accelerators  for  high-energy  physics  research  are  to  be  realized.  To  this  end,  the  primary 
thrust  of  this  research  is  the  development  of  a  model  for  a  stacked  Blumlein  microwave  source  used  in  such 
a  radar  system. 

Model  Development 

In  1948,  A.  D.  Blumlein  patented  a  device2  capable  of  both  pulse  formation  and  propagation.  It 
consists  of  two  transmission  lines  that  share  a  common,  charged  conductor.  When  connected  to  a  matched 
load  impedance,  ZL,  the  measured  output  is,  ideally,  a  rectangular  pulse  with  an  amplitude  that  is  equal  to 
the  full  charging  voltage,  V0.  The  width  of  this  pulse  is  roughly  twice  the  one-way  transit  time  of  a 
constituent  transmission  line.  The  remarkable  achievement  of  this  invention  is  that  it  overcomes  the 
obstacle  of  a  regular  transmission  line  for  which  only  one  half  of  the  charging  voltage,  at  best,  is  observable 
across  a  matched  load.  Because  the  transmission  lines  are  effectively  charged  in  parallel  and  later 
discharged  in  series,  a  stack  of  Blumleins  (Figure  1)  can  be  constructed  whose  output  approaches  the 
charging  voltage  multiplied  by  the  number  of  Blumleins  comprising  the  stack,  thereby  producing  the 
additional  benefit  of  voltage  multiplication. 

The  charging  and  discharging  of  Blumleins  is  essentially  a  time  domain  problem.  While  the 
propagation  of  a  pulse  on  a  lossy  line  is  best  dealt  with  in  the  frequency  domain,  the  transmission  lines 
considered  in  this  work  are  lossless.  This  means  that  frequency  dependent  dispersion  is  not  a  factor  as  long 
as  only  a  single  dielectric  is  present.  The  computational  analysis,  therefore,  utilizes  space-time  diagrams,3 
otherwise  known  as  bounce  diagrams,  for  a  purely  time-domain  treatment  of  stacked  Blumlein  pulse 
generators.  However,  the  method  of  moments  is  first  employed  to  calculate  the  system  parameters,  such  as 
the  matrices  for  the  capacitance  per  unit  length  and  the  characteristic  impedance. 
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Figure  1:  A  stack  of  two  Blumleins  in  the  (A)  shared  and  (B)  separate  conductor  configurations.  In  both 
instances,  the  conducting  plates  are  all  separated  by  a  distance  d.  For  configuration  in  (B),  however,  the 

spacing  between  Blumleins  is  variable. 


In  this  research,  the  moment  method  was  applied  in  two  dimensions  by  employing  the  following 


equation 
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where  0,  is  the  known  potential  at  the  center  of  the  rth  subsection.  The  first  term  on  the  right  hand  side  of 
the  equation  represents  the  potential  of  the  rth  subsectional  element  due  to  itself.  The  second  term  reflects 
the  influence  of  all  other  elements.  In  this  case,  w  is  the  width  of  each  subsection.  It  is  constant  for  all  M 
elements  of  the  structure.  In  other  words,  w,  =  wj  =  w  =  W/m.  The  variable  X,  is  the  linear  charge  density 
for  each  element  to  be  determined  in  the  computation.  Although  it  is  evenly  distributed  over  each  element, 
the  charge  density  is  treated  as  a  point  charge  located  at  the  center  of  the  subsection.  The  distance  between 
elements  is  r,  -  r}.  Thus  the  final  form  of  equation  above  is 
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Figure  2:  Cross-sectional  view  of  a  stack  of  N  Blumleins.  Each  conducting  plate  is  divided  into  m  equal 
segments,  creating  a  total  of  M  elements  for  the  entire  structure. 


There  are  M  of  these  equations  to  solve  simultaneously  for  the  M  charge  densities.  This  is  easily 
accomplished  with  a  computer  software  package  designed  to  handle  matrix  operations,  such  as  MATLAB. 
The  matrix  equation  used  for  this  research  is 

M-icra 

where  [O]  is  an  M  x  1  matrix  of  the  potential  on  each  subsectional  element.  [C]'1  is  the  inverse  of  the 
capacitance-per-unit-length  matrix.  The  inverse  capacitance  matrix  can  be  generalized  as  follows 

(cr—t&i 

£ 

where  [fg]  is  the  “per-element”  matrix  of  geometric  factors.5  This  matrix  is  computed  instead  of  [C]'1 
because  it  can  be  used  to  calculate  matrices  for  quantities  such  as  inductance,  impedance  and  admittance. 
In  order  to  proceed  with  the  development  of  the  model,  however,  the  “per-element”  matrix  of  geometric 
factors  has  to  be  reduced  to  a  matrix  with  “per-plate”  entries,  with  one  of  those  plates  as  a  reference.  This  is 
required  because  the  electric  potential  between  plates  is  measured  with  respect  to  some  point  of  reference. 
Without  this  reduction,  we  would  therefore  be  unable  to  properly  keep  track  of  the  transient  waves  (and 
their  associated  voltages)  established  and  propagated  by  a  stack  of  Blumleins. 

The  process  of  reducing  the  “per-element”  matrix  of  geometric  factors  to  a  “per-plate”  one  is  not 
accomplished  by  simply  summing  all  [fg]  matrix  elements  for  each  plate.  When  calculating  the  matrix  of 
geometric  factors  per  plate,  [FG],  we  must  first  invert  [fg],  then  perform  the  summation  of  [gf]  matrix 
entries  for  each  plate.  This  produces  [GF],  which  is  the  n  x  n  inverted  “per-plate”  matrix  of  geometric 
factors  having  matrix  elements 
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Inverting  the  [GF]  matrix  finally  yields  the  “per-plate”  matrix  of  geometric  factors.  Now,  this  needs  to  be 
reduced  yet  further  so  that  the  bottom-most  plate  becomes  the  point  of  reference  for  the  system.  This  is 
accomplished  with  the  following  equation 


Fgy  =  FGmj+x  -  FGi+u  -  FGi  J+l  +  FGU  i  =  1 ...«  - 1  j  =  1 ...»  - 1 , 

where  Fgy  are  the  elements  for  the  reduced  “per-plate”  geometric  factor  matrix,  [Fg]. 

Now  that  these  fundamental  matrices  have  been  determined,  it  is  straightforward  to  compute  the 
capacitance  and  characteristic  impedance  of  the  stacked  Blumlein  structure.  The  per-unit-length 
characteristic  impedance  matrix  is  given  by 


while  the  other  circuit  parameter,  the  per-unit-length  capacitance,  is  related  to  [Fg]  by 

[< cMj%: r. 


These  matrices  alone  do  not  give  all  of  the  information  needed  to  track  the  propagation  of  TEM 
waves.  They  do,  however,  prove  necessary  in  the  determination  of  the  scattering  matrix  for  each 
termination  of  the  stacked  Blumlein.  At  the  switching  end,  we  find 

(]z n]-[ap 
K[zn\*[zc\)’ 

where  [ZT1]  is  the  termination  impedance  matrix  at  the  switch  end  of  the  Blumlein.  At  the  opposite 
termination  we  find 


N= 


^[Z7’2]-[Zc]> 

v[ZT2]+[ZcL 


[ZT2]  is  the  termination  impedance  matrbc  for  the  load  end  of  the  Blumlein.  Consequently,  a  properly 
matched  load  must  be  determined  that  is  a  scalar  value,  not  a  matrix. 

As  seen  in  Figure  1,  the  load  is  placed  between  the  outermost  plates.  Ideally,  a  load  whose 
impedance  matches  that  of  the  Blumlein  stack  will  have  an  optimized  output.  As  Figure  3(A)  indicates,  the 
load  “sees”  the  series  summation  of  the  four  adjacent-plate  impedances  if  there  is  no  field  coupling  between 
non-adjacent  plates.  Complications  arise  because  that  field  coupling  does  indeed  exist,  as  evidenced  in 
Figure  3(B). 

The  stacked  Blumlein  model  computes  and  suggests  two  reasonably  “matched”  loads.  The  first  is 
simply  the  series  summation  of  each  transmission  line’s  characteristic  impedance  (Figure  3(A)) 
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A)  B) 


Figure  3:  (A)  The  generally  accepted  view  of  the  characteristic  impedances  involved  in  obtaining  a 
“matched”  load  of  a  stacked  Blumlein.  (B)  The  matrix  formulation  of  this  problem  reveals  all  of  the 
“other”  characteristic  impedances  present  for  this  structure,  complicating  the  computation  of  a  true  “match.” 
(The  placement  of  impedances  does  not  reflect  their  actual  position.) 


The  second  includes  coupling  contributions.  For  a  stack  of  two  Blumleins  that  have  no  spacing  between 
them  (Figure  3(B)),  this  amounts  to 
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Since  the  load  is  connected  solely  between  the  outermost  conductors,  we  need  only  concern  ourselves  with 
single  or  equivalent  series  impedances  that  mirror  that  connection.  * 

Having  determined  all  of  the  pertinent  matrices,  all  that  remains  is  the  utilization  of  space- 
time  diagrams  to  track  the  propagation  of  the  transient  TEM  waves.  This  is  accomplished  with  voltage 
vectors  and  matrices.  As  Figure  4  indicates,  all  potential  differences  are  measured  with  respect  to  the 
bottom-most  plate.  In  this  case,  the  electrostatic  potentials  are  shown.  If  the  stack  of  Blumleins  shown  in 
Figure  4  is  charged  to  a  voltage  V0,  the  resulting  static  potential  vector  is 

v; 
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Figure  4:  The  “illustrated”  components  of  the  static  potential  vector.  “Vector”  simply  refers  to  a  single 

column  matrix. 


and  the  static  voltage  matrix  is 

V0  -V0  0  -V0~ 

-V0  0  V0  0 

0  v0  V0  -V0  • 

-v0  o  -V0  0  _ 

Once  the  switches  are  closed,  the  first  transient  waves  are  established, 
determined  by 


M= 


~\Zc]  ) 

Uzc]+[zn]J 


N- 


Their  voltage  component  is 


When  the  waves  reach  the  load  termination,  they  undergo  another  reflection  governed  by  the  scattering 
matrix  [SLJ.  Thus,  the  second  set  of  traveling  waves  have  voltage  components  evaluated  by 

[Vt2]  =  [SL\Vti\. 


(Although  not  shown  here,  the  total  voltage  vectors  and  matrices  for  each  termination  are  also  computed.) 
During  the  third  transit  time,  the  second  set  of  transient  waves  reaches  the  z  =  0  termination,  where  the 
scattering  matrix  is  [SO]  and  forms  the  third  set  of  traveling  waves  whose  voltage  components  are 

[Vt3]=[S0lVt2], 

The  process  is  repeated  for  ten  transit  times,  which  is  usually  long  enough  to  establish  the  dissipation  of 
energy  for  well  matched  loads. 
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Results 


The  simulation  of  pulse  formation  and  propagation  in  stacked  Blumleins  is  performed  by  a 
program  written  with  MATLAB  (Version  4.0  for  Windows).  This  analysis  is  very  versatile  despite  being 
limited  to  a  single  homogeneous  dielectric  and  only  one  basic  conducting  plate  geometry.  The  flexibility 
arises  from  the  number  and  type  of  variables  that  are  required  as  input.  First  and  foremost,  the  user  must 
enter  the  number  of  Blumleins  being  stacked.  The  code  is  capable  of  handling  a  single  Blumlein,  as  well  as 
stacks  of  two  or  more.  If  the  user  wishes  to  investigate  multiple  Blumleins,  they  are  then  prompted  for  the 
vertical  distance  between  the  Blumleins,  where  an  entry  of  zero  corresponds  to  shared  conductors,  as  seen 
in  Figure  1. 

The  user  must  enter  the  length  and  width  of  the  constituent  conducting  plates  of  the  Blumleins,  as 
well  as  the  distance  that  separates  them.  The  thickness  of  a  plate  is  considered  infinitesimally  thin.  Each 
plate  in  the  stack  is  assumed  to  have  identical  dimensions  and  equal  vertical  spacing.  The  volume  between 
and  around  the  conductors  is  filled  with  dielectric.  The  relative  permittivity  is,  therefore,  the  final  bit  of 
structural  information  requested  as  input. 

The  accuracy  of  the  simulation  depends  on  the  number  of  elements  into  which  a  cross-sectional 
slice  of  the  stack  is  divided.  For  instance,  reasonable  results  can  be  obtained  if  a  conductor  with  a  width  of 
2.54  cm  is  sectioned  into  15  elements.  The  code,  in  that  case,  takes  a  noticeably  shorter  amount  of  time  to 
run  than  if  50  elements  were  chosen.  Of  course,  better  results  are  generated  as  the  number  of  elements 
increases,  but  this  is  at  the  expense  of  increased  computational  time.  The  choice  of  acceptable  accuracy  is 
left  to  the  discretion  of  the  user.  For  all  plots  displayed  here,  however,  a  2.54-cm  wide  Blumlein  is  always 
divided  into  50  subsections. 

The  user  must  also  input  the  operational  parameters  of  the  stacked  Blumlein  system.  These  include 
the  charging  voltage  and  the  load  impedance.  As  mentioned  earlier,  two  different  load  impedance  values 
are  calculated  and  displayed  on  the  screen.  The  user  can  choose  an  impedance  that  optimizes  either  the 
voltage  multiplication  factor  or  the  power  of  the  main  output  pulse,  or  opt  for  another  value  entirely. 

Due  to  errors  that  will  arise  (division  by  zero/multiplication  by  infinity),  every  open  and  short  are 
assigned  finite,  non-zero  values.  Since  this  simulation  is  wont  to  treat  a  closed  switch  as  a  short,  the  user,  is 
presented  with  the  opportunity  of  determining  the  resistance  of  the  switching  element.  Although  this 
treatment  is  not  entirely  accurate,  it  is  more  than  adequate  for  this  research,  for  the  creation  of  a  thorough 
model  of  a  semiconducting  photocondutive  switch  or  a  thyratron  is  beyond  the  scope  of  this  work. 

By  requiring  all  of  the  variables  mentioned  above,  the  model  is  more  complete  than  if  the  user 
were  only  to  input  the  number  of  Blumleins,  the  characteristic  impedance  of  the  Blumleins,  and  the 
charging  voltage  for  the  system.  Even  though  a  single  characteristic  impedance  is  measured  in  the  lab, 
theory  dictates  that  it  is  produced  by  a  combination  of  characteristic  impedance  measurements.  Therefore, 
this  computer  model  of  stacked  Blumleins  incorporates  all  of  this  information  to  determine  the  capacitance 
and  characteristic  impedance  of  the  structure.  These  matrices  are  then  used  to  calculate  the  initial  energy 
stored  by  the  Blumlein  and  the  “matched”  load  impedances,  respectively.  The  initial  energy  is,  in  turn, 
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utilized  in  the  computation  of  the  energy  efficiency  of  the  stacked  Blumlein.  The  load  impedance  is  fed  into 
the  portion  of  the  program  responsible  for  generating  the  scattering  matrix  for  the  z  =  /  termination,  and  also 
into  the  equations  used  to  compute  the  final  energy  seen  across  the  load. 

Because  of  this,  the  matrices  for  characteristic  impedance,  capacitance,  and  scattering  at  both  ends 
of  the  stacked  Blumlein  are  displayed  on  the  screen.  It  is  trivial  to  modify  the  code  so  that  it  displays  any 
other  matrix  or  vector  encountered  during  each  run  of  the  program.  This  might  include  the  voltage  vector 
for  any  of  the  transient  waves  or  the  matrix  of  total  voltage  between  all  the  plates  in  the  stack  for  either 
termination  at  any  desired  time.  All  of  the  pertinent  information  about  the  structure  and  operational 
parameters  of  the  Blumleins  is  available  on  the  generated  plots  of  voltage  against  transit  time,  which  can  be 
printed  to  the  default  printer.  Perhaps  the  information  of  most  interest  that  appears  on  these  graphs  is  the 
energy  efficiency  and  the  voltage  multiplication  factor.  The  efficiency  is  determined  by  dividing  the  energy 
in  the  output  pulse  by  the  total  energy  initially  stored  by  the  stacked  Blumleins,  or 
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if,  for  example,  there  are  two  Blumleins,  and 

v2 

^  ou?, max  ~ 

Efmal  =  zZ  2t 

The  energy  initially  stored  in  the  structure  is  calculated  by  adding  the  energy  stored  between  all  plates. 
Obviously,  some  sets  of  plates  store  nothing  because  the  electrostatic  potential  between  them  is  zero.  The 
output  energy  is  comparatively  easy  to  compute  as  it  involves  only  the  maximum  voltage  of  the  main  output 
pulse,  the  temporal  width  of  that  pulse,  and  the  load  impedance. 

The  voltage  multiplication  factor  is  calculated  by  dividing  the  maximum  output  voltage  of  the  main 
pulse  by  the  charging  voltage. 


Ideally,  the  voltage  multiplication  factor  should  be  equal  to  the  number  of  Blumleins  being  stacked.  As  the 
following  simulation  results  indicate,  this  hardly  ever  occurs.  What  is  noticeable,  though,  are  the  trends  that 
appear  with  changes  to  input  variables  such  as  the  number  of  Blumleins  and  the  choice  of  load  impedance. 
As  can  be  deduced  from  the  plot  below,  there  is  a  definite  correlation  between  efficiency  and  the  number  of 
Blumleins.  It  is  easy  to  conclude  that  as  the  number  of  Blumleins  increases,  the  efficiency  suffers  a  marked 
decrease.  The  graph  in  Figure  5  shows  this  trend  as  the  ratio  of  plate  separation,  d,  to  plate  width,  w,  varies. 
We  might  gather  from  the  plot  that  a  more  efficient  machine  can  be  built  by  minimizing  the  geometric  ratio. 
In  reality,  though,  there  is  a  limit  as  to  what  is  physically  plausible.  The  graph  also  suggests  that  the  number 
of  Blumleins  in  a  vertical  stack  should  be  kept  below  four  or  possibly  five  due  to  the  fact  that  beyond  this 
limit,  the  efficiency  drops  below  what  most  would  consider  an  acceptable  level. 
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Figure  5:  Variation  of  efficiency  with  the  number  of  Blumleins  in  a  stack  and  the  geometric  ratio,  d/w . 
The  following  variables  are  constant:  dielectric  constant  (transformer  oil),  Blumlein  separation  (none),  and 

load  impedance  (efficiency  maximized). 

The  trends  of  the  graph  in  Figure  5  are  hardly  surprising.  As  the  geometric  ratio  increases,  the 
capacitance  decreases  and  the  characteristic  impedance  increases.  Clearly,  the  amount  of  energy  initially 
stored  will  be  affected  (diminished)  by  the  change  in  capacitance.  The  variation  in  the  characteristic 
impedance  translates  to  an  alteration  of  the  computed  load  impedance  suggestions.  As  the  characteristic 
impedance  increases,  so  does  the  load  impedance  since  it  is  calculated  with  the  elements  of  the 
aforementioned  matrix.  Recall,  too,  that  the  load  impedance  is  used  to  determine  the  termination 
impedance  matrix,  which  is  in  turn  involved  in  the  evaluation  of  the  scattering  matrix  for  the  load  end  of  the 
Blumleins.  Additionally,  the  characteristic  impedance  matrix  is  itself  a  component  of  the  scattering  matrix 
calculation.  Thus,  a  change  in  the  elements  of  the  scattering  matrix  will  evoke  a  response  in  the  vectors  and 
matrices  of  the  voltage  component  of  the  waves  formed  and  propagated  by  the  Blumleins.  As  the  maximum 
voltage  across  the  load  diminishes,  so  too  will  the  energy  and  efficiency.  Clearly,  then,  increasing  the 
geometric  ratio  only  serves  to  decrease  the  efficiency  of  a  stack  of  Blumleins. 

As  the  number  of  Blumleins  increases,  field  fringing  and  coupling  occurs  between  more  and  more 
conductors.  The  strength  of  field  coupling  is  related  to  the  characteristic  impedance.  The  strongest 
coupling  (smallest  value  of  Zc)  occurs  between  adjacent  plates.  Hence  the  majority  of  energy  is  found  in 
the  waves  propagating  between  these  plates.  Even  though  the  coupling  associated  with  non-adjacent  plates 
is  weaker  than  that  for  adjacent  conductors,  the  overall  effect  is  the  reallocation  of  some  energy  that  would 
have  otherwise  been  propagated  between  the  adjacent  conductors.  We  can  infer,  then,  that  the  voltage 
associated  with  the  propagating  waves  likewise  diminishes.  Thus,  the  efficiency  of  the  device  is  lowered. 
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Figure  6:  Variation  of  the  voltage  multiplication  factor  with  the  number  of  Blumleins  in  a  stack  and  the 
geometric  ratio,  d/w.  The  following  variables  are  constant:  dielectric  constant  (transformer  oil),  Blumlein 
separation  (none),  and  load  impedance  (efficiency  maximized). 


Parallel  conclusions  can  be  drawn  upon  investigation  of  the  plot  of  the  voltage  multiplication 
factor  versus  the  number  of  Blumleins  (Figure  6).  Because  the  maximum  output  voltage  decreases  with 
increasing  d/w ,  the  voltage  multiplication  factor  will  also  suffer  the  same  fate.  Once  again,  the  model  seems 
to  suggest  that  stacks  of  no  more  than  four  Blumleins  are  worth  constructing  in  this  particular  situation. 

The  impact  of  conductor  spacing  is  clear.  Now  consider  the  effect  of  the  spacing  between 
Blumleins,  as  depicted  in  Figure  7.  We  would  expect  the  efficiency  to  fall  off  as  the  Blumlein  separation 
increases  for  the  same  reasons  discussed  when  the  conductor  separation  was  increased.  Basically,  the  losses 
arising  from  the  coupling  of  fringing  fields  become  more  pronounced  with  larger  Blumlein  separation 
because  an  additional  wave  is  propagated  between  the  top-most  conductor  of  the  lower  Blumlein  and  the 
bottom-most  conductor  of  the  upper  Blumlein.  The  wave  is  obviously  eliminated  when  these  two 
conductors  are  replaced  by  a  single,  shared  plate.  No  energy  is  initially  stored  between  these  two 
conducting  plates  (they  are  at  the  same  potential),  so  increasing  the  separation  between  them  does  not  affect 
the  amount  of  energy  initially  stored  in  the  stack  of  Blumleins.  The  energy  of  the  wave  that  travels  between 
the  conductors  after  switch  closure  must  therefore  detract  from  the  energy  available  to  the  waves  propagated 
by  other  adjacent  conductors,  thereby  lowering  the  efficiency  of  the  device. 

The  marked  difference  in  efficiencies  does  not  appear  when  the  voltage  multiplication  factor  is 
plotted  against  the  number  of  Blumleins  for  a  variation  in  Blumlein  separation  (Figure  8).  Instead, 
increasing  the  separation  between  the  Blumleins  apparently  increases  the  voltage  multiplication  factor.  This 
comes  about  as  a  result  of  the  extra  wave  propagated  between  each  pair  of  Blumleins.  Its  voltage  adds  to  all 
of  the  others,  increasing  the  output  despite  the  fact  that  all  of  the  other  waves  have  a  slightly  lower  voltage 
as  a  result  of  the  presence  of  that  wave. 
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Figure  7:  Variation  of  the  efficiency  with  the  number  of  Blumleins  in  a  stack  and  the  Blumlein  separation, 
blsep.  The  following  variables  are  constant:  dielectric  constant  (transformer  oil),  geometric  ratio 
( d/w  =  0.3),  and  load  impedance  (efficiency  maximized). 
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Figure  8:  Variation  of  the  voltage  multiplication  factor  with  the  number  of  Blumleins  in  a  stack  and  the 
Blumlein  separation,  blsep.  The  following  variables  are  constant:  dielectric  constant  (transformer  oil), 
geometric  ratio  ( d/w  =  0.3),  and  load  impedance  (efficiency  maximized). 

The  next  two  plots  (Figures  9  and  10)  nicely  illustrate  the  effect  that  different  load  impedances 
have  on  the  efficiency  and  voltage  multiplication  factor.  Earlier,  we  saw  the  formulation  of  a  “matched” 
load  by  omitting  or  including  the  characteristic  impedances  of  some  of  the  non-adjacent  plates.  It  was 
suggested  that  a  truer  match  might  be  made  by  including  more  characteristic  impedances.  In  other  words, 
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the  latter  formulation  was  more  accurate  because  it  allowed  a  larger  portion  of  the  wave  to  be  transmitted, 
which  meant  that  there  were  fewer  reflections  and  that  the  stored  energy  was  dissipated  more  quickly. 

Since  the  load  impedance  is  determined  from  the  characteristic  impedance  matrix,  the  inclusion  of 
more  characteristic  impedances  in  the  load  impedance  computation  results  in  a  smaller  calculated  load 
impedance.  Clearly,  this  also  affects  the  scattering  matrix  for  the  load  end  of  the  Blumlein,  thereby 
decreasing  the  maximum  output  voltage.  Since  the  change  in  voltage  is  less  pronounced  than  the  change  in 
load  impedance,  a  greater  efficiency  is  achieved.  Thus,  this  method  of  determining  a  “match”  is  referred  to 
as  optimizing  for  higher  output  power  or  greater  energy  efficiency. 

The  omission  of  the  characteristic  impedances  of  non-adjacent  plates,  on  the  other  hand,  results  in 
an  increased  load  impedance  because  nothing  is  in  parallel  with  the  series  equivalent  of  the  adjacent  plate 
impedances.  This,  then,  decreases  the  efficiency  because  the  output  energy  is  divided  by  a  larger  number, 
despite  the  small  increase  in  the  maximum  output  voltage.  It  does  clearly  increase  the  voltage 
multiplication  factor  (Figure  10)  because  of  that  increase  in  output  voltage.  Hence,  this  “match”  is  denoted 
as  optimization  for  a  higher  voltage  multiplication  factor.  Obviously,  if  we  were  designing  a  stacked 
Blumlein  pulse  generator,  we  would  have  to  decide  which  is  more  important  for  our  intended  application,  a 
device  which  has  better  voltage  multiplication  capabilities  or  one  that  produces  higher  power  pulses. 


Number  of  Blumleins 


Figure  9:  Variation  of  the  efficiency  with  the  number  of  Blumleins  in  a  stack  and  the  load  impedance,  ZL 
(energy  efficiency  or  voltage  multiplication  number  optimized).  The  following  variables  are  constant: 
dielectric  constant  (transformer  oil),  geometric  ratio  (d/w  =  0.3),  and  Blumlein  separation  (none). 
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Number  of  Blumleins 

Figure  10:  Variation  of  the  voltage  multiplication  factor  with  the  number  of  Blumleins  in  a  stack  and  the 
load  impedance,  ZL  (energy  efficiency  or  voltage  multiplication  factor  optimized.  The  following  variables 
are  constant:  dielectric  constant  (transformer  oil),  geometric  ratio  (d/w  =  0.3),  and  no  Blumlein  separation. 

Conclusions 

(Due  to  the  limited  amount  of  space  available,  this  section  will  be  kept  unusually  brief.  For  a  more  in-depth  discussion,  please  refer 
to  my  dissertation.) 

The  model  of  stacked  Blumleins  developed  for  this  research  was  verified  using  three  different 
procedures.  The  matrices  for  the  capacitance,  characteristic  impedance,  and  scattering  were  checked 
against  a  frequency-domain  software  package  from  France  called  CRIPTE.6  Excellent  agreement  was 
exhibited  between  the  two  models.  The  efficiency  was  compared  to  a  QBASIC  analysis  written  by  lab 
scientist  Jim  O’Loughlin7  that  approximated  a  stack  of  Blumleins  as  a  network  of  resistors.  Due  to  major 
differences  in  treatment  of  this  problem  (circuit  theory  versus  electromagnetic  wave  theory),  agreement  was 
not  as  well  defined,  although  both  methods  concur  that  constructing  stacks  of  more  than  two  to  four 
Blumleins  would  be  unwise  because  of  operation  inefficiency.  Finally,  experimental  results  for  a  stack  of 
two  Blumleins  were  used.  Again,  the  model  agreed  with  reality  as  well  as  could  be  expected  considering 
the  approximations  made  during  code  development. 

Obviously,  more  work  can  still  be  done.  Future  extensions  of  this  research  might  include  the 
introduction  of  multiple  dielectrics,  conducting  plates  of  variable  width,  or  conductor  thickness.  A  more 
rigorous  treatment  of  the  terminations  should  also  be  attempted.  For  example,  non-simultaneous  switch 
commutation  might  be  investigated,  or  perhaps  the  effects  of  attaching  a  TEM  horn  antenna  might  be 
explored.  This  naturally  leads  to  questions  pertaining  to  the  determination  of  a  better  matched  load. 
Finally,  different  parallel  plate  Blumlein  arrangements  should  be  considered,  say  side-by-side  instead  of 
vertical  stacking.  Clearly  improvements  can  be  made,  but  the  model,  as  it  currently  stands,  does  provide  a 
reasonable  starting  point  for  a  researcher  designing  a  stacked  Blumlein  pulse  generator. 
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Abstract 

Using  the  results  derived  by  Akira  Matsumoto  and  Kenichi  Iwamoto 
published  in  two  papers12  studying  Morse  potentials,  we  attempted  to  calculate 
transition  probabilities  between  vibrational  levels.  We  used  two  methods  for 
calculating  the  probabilities,  the  first  method  involved  numerically  evaluating  the 
analytic  eigenfunctions  and  the  second  method  was  to  evaluate  the  analytic 
solution  to  the  electric  dipole  matrix  elements,  both  of  which  were  provided  by 
Matsumoto  and  Iwamoto’s  papers. 
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Introduction 


Much  research  has  been  done  on  electric  dipole  matrix  elements  and 
vibrational  transition  probabilities  and  much  data  has  been  provided  to 
approximate  these  values.  However,  there  is  still  more  to  be  learned  at  high 
vibrational  levels  and  there  is  still  a  need  for  better,  more  accurate  models  of  the 
potentials  of  diatomic  molecules.  We  attempted  to  get  numerical  values  for  the 
current  analytical  solutions  to  then  use  for  analysis  of  spectra  obtained  in  the 
laboratory. 


Methodology 

For  the  first  approach,  we  wrote  computer  programs  in  Series  Processing 
Language  for  use  with  Dadisp  software.  The  code  was  written  in  a  style  similar 
to  C  and  the  purpose  of  the  programs  was  to  numerically  calculate  the  solutions 
to  the  analytic  eigenfunctions.  The  eigenfunctions  contain  associated  LaGuerFe 
polynomials,  so  we  had  to  write  code  to  evaluate  these  polynomials  using  the 
recurrence  function  for  these  polynomials.  Once  the  numeric  results  for  the 
eigenfunctions  were  obtained,  we  plotted  the  Morse  potential  (using  the  equation 
in  the  two  referenced  papers)  and  overlaid  the  eigenfunctions  (which  were 
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normalized)  with  the  potential  and  saw  that  they  were  well  behaved.  The  next 
step  in  this  project  is  to  numerically  calculate  the  transition  probability  between 
two  levels  by  integrating  the  product  of  the  two  eigenfunctions  and  the  radius 
(the  operator).  We  expect  that  these  numbers  will  match  the  numbers  obtained 
in  our  other  approach  to  the  problem. 

Our  other  approach  was  to  take  the  analytic  solution  to  the  electric  dipole 
matrix  elements  derived  by  Iwamoto  and  Matsumoto1’2  and  numerically  evaluate 
it  for  molecular  nitrogen  using  experimental  spectroscopic  values  for  We,  WeXe 
and  re  taken  from  the  Journal  of  Physical  and  Chemical  Reference  Data3.  We 
confirmed  that  the  programs  written  to  carry  out  the  calculations  were  working 
correctly  by  first  reproducing  the  data  for  CO  reported  by  Matsumoto  and 
Iwamoto1,2  before  we  then  changed  the  parameters  for  our  own  experimental 
interests  (N2).  We  were  then  able  to  plot  this  data  as  relative  transition 
probability  versus  the  vibration  level  and  see  the  trends  in  probability  of 
transitions  as  you  go  up  in  vibrational  level.  The  next  step  in  this  project  is  to 
see  if  these  numbers  match  the  probabilities  calculated  using  the  original 
eigenfunctions  and  to  apply  this  information  to  experimental  data.  By 
theoretically  modeling  the  transition  probabilities,  we  can  analyze  experimental 
spectra  by  comparison  of  theoretical  predictions  to  experimental  results. 
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Results 


At  this  point,  it  would  be  preliminary  to  report  any  official  results  because 
the  two  phases  of  the  project  have  not  yet  been  compared.  Once  it  can  be 
confirmed  that  both  approaches  to  solving  the  problem  give  logical  results,  then 
these  values  can  be  reported  and  utilized  in  analysis  of  experimental  spectra. 

However  it  can  be  stated  that  at  this  point,  all  codes  are  running  properly 
and  appear  to  be  giving  correct  numerical  results.  This  is  evidenced  by 
overplotting  the  eigenfunctions  with  the  Morse  potential  and  seeing  that  they  are 
well  behaved  as  well  as  seeing  reasonable  looking  plots  of  relative  transition 
probability  versus  v  level.  This  project  is  still  underway  and  hopefully  there  will 
be  some  quantitative  results  to  report  in  the  near  future. 


Discussion 


One  suggestion  for  improvement  to  this  experiment  that  may  be 
implemented  before  its  conclusion  is  to  utilize  a  more  powerful  and  scientifically 
capable  programming  language  to  carry  out  the  calculations  for  actual  data 
analysis.  There  still  exists  some  concern  that  round  off  error  and  the  limitations 
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on  variable  size  may  cause  deviation  from  actual  values.  This  is  an  area  still 
under  investigation.  The  programs  involve  very  large  numbers  and  exponents 
as  well  as  utilize  complex  mathematical  functions  such  as  the  gamma  function 
and  the  recurrence  functions  for  associated  LaGuerre  polynomials.  Therefore 
the  values  of  some  variables  within  the  program  often  approach  the 
technological  limits  of  the  system  on  which  they  are  being  run. 

However  our  results  at  this  point  in  the  project  are  reasonably  similar  to 
those  found  by  other  researchers.  This  is  evidenced  by  the  ability  of  our 
programs  to  reproduce  the  electric  dipole  matrix  element  values  for  CO  reported 
by  Iwamoto  and  Matsumoto1,2. 


Conclusion 


At  this  point  in  the  project,  we  feel  that  we  have  developed  accurate  and 
functional  programs  to  calculate  Morse  potentials,  vibrational  transition 
probabilities  and  vibrational  electric  dipole  matrix  elements.  We  intend  to 
finalize  these  programs,  confirm  the  accuracy  of  their  results  and  utilize  this 
information  in  the  analysis  of  spectroscopic  data  taken  in  the  laboratory. 
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CONDUCTING  FLUID  REAL-TIME  2-D  ELECTRONIC  INTERPOLATOR  AND  SPATIAL 
FILTER  FOR  WAVEFRONT  SENSOR-TO-CORRECTOR  INTERFACING. 


Tyrone  Ospino-Marthe 

Graduate  Student 
Department  of  Physics 
University  of  Puerto  Rico/Mayaguez  Campus 


INTRODUCTION 

The  situation  often  exist  in  adaptive  optics  that  the  geometry  of  wavefront  sensing  is  different 
than  that  of  wavefront  correcting.  By  virtue  of  Nyquist  Theorem,  more  wavefront  sensing 
elements  are  usually  required  than  for  wavefront  correcting,  and  the  geometrical  arrangement 
may  differ.  In  our  case,  the  wavefront  sensor  has  rectangular  geometry  while  the  correcting 
element  has  hexagonal.  Another  situation  arises  when  the  reconstructed  wavefront  solution  is 
noisy  and  must  be  spatially  filtered  to  the  lower  order  modes.  We  propose  a  solution  using  the 
Fluidic  Interpolator  Modulefhenceforth  called  FILM). 

An  aqueous  solution  of  CuS04  provides  for  a  resistive  medium  which  allows  for  freedom  of 
movement.  Electronically  driven,  etched  copper  circuit  board  traces  and  pads  with  the  input 
geometry  provide  for  spatially  defined  potentials  in  the  medium.  A  receiving  plate  etched  with 
the  field  distribution  and  resulting  fluid  potentials  interpolates  and  smoothes  the  output  plate 
data.  When  the  output  plate  is  “close”  to  the  input  plate  i.e.  the  separation  is  much  less  than  the 
separation  between  input  plate  elements,  the  output  plate  senses  essentially  bilinear 
interpolation.  Separating  the  plates  tends  toward  averaging  of  the  fields,  and  hence  2-D 
smoothing  and  interpolation  of  the  output  data. 
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CONDUCTING  FLUID  REAL-TIME  2-D  ELECTRONIC  INTERPOLATOR  AND  SPATIAL 
FILTER  FOR  WAVEFRONT  SENSOR-TO-CORRECTOR  INTERFACING. 


Tyrone  Ospino-Marthe 


Methodology 

The  procedures  for  solving  these  problems  are  the  following: 

a.  I  need  to  calculate  in  every  point  in  the  square  16x16  the  capacitance,  next  calculate  the 
contribution  to  every  point  to  the  square  on  the  hexagonal,  and  to  know  the  values  for  the 
capacitance  in  every  127  points  in  the  hexagon. 

b.  Assuming  and  entry  to  voltages  in  the  wavefront,  calculate  the  total  voltages  in  the  square. 

C.  Calculate  next  the  Charge  in  every  point  to  the  square  128  X  128. 

d.  Calculating  the  first  layer,  make  it  16  times  and  calculating  the  16  layers. 

e.  Addition  noise  to  the  principal  function  and  make  the  process  other  time. 

f.  After  to  calculate  the  first  layer,  calculate  again  the  16  layers. 

g.  Make  the  electronic  part  constructing  two  parallels  boards  ,  one  square  128x128  and  other 
hexagesimal  127  under  a  “place”  fill  with  CUS04+5H20,  with  the  three(3)  axis  movements. 


Results 


Here  we  have  to  solution  many  problems,  one  to  the  principal  restriction  are  that  I  need  to  work 
in  three-dimensions,  and  the  formulas  and  bibliography's  don’t  mentioned  this  cases  in  his 
examples.  First,  need  us  to  solve  and  have  clear  the  physics  process  here. 


> 

We  idealize  an  square  to  16x16  and  apply  this  formulas  for  the  corners,  sides  and  interior. 


CORNERS 

Qii  =  qo(2+1/V2) 

SIDES 

q„=q0(3+2/V2) 

INTERIOR 

q„=q0(4+4/V2) 

must  be  remember  ,  in  this  point  the  calculate  is  to  the  relative  capacitance  not  the  “real” 
capacitance,  with  this  values  we  making  a  CAPACITANCE  MAP  (  in  all  this  process,  the 
simulations,  I  will  make  it  with  matlab’s  programs). 
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Next  to  this  calculus,  find  the  contribution  for  every  point  to  the  square  over  the  hexagon, 
calculate  in  every  127  point  to  the  hexagon  with  127  formulas,  seeing  geometrically  the 
contribution  with  a  superposition  to  the  hexagon  over  the  square  and  neglected  the  comers  . 
Maybe  one  way  to  make  this  calculus  are  using  numerical  methods,  maybe  the  relaxation 
method,  images  method,  or  other ,  but  we  can’t  used  for  geometrical  problems. 

Now  knowing  the  contribution  into  the  boards  ,  we  continuos  with  a  program  for  to  make  a  net 
with  only  voltages  in  specifics  points  in  a  matrix  to  128x128,  the  positions  are:  (row  .column) 


1,1 

1,9 

1,18 

1,27 

... 

1,126 

1,128 

9,1 

9,9 

9,18 

9,127 

9,128 

18,1 

18,9 

18,18 

18,126 

18,128 

128,1 

128,9 

128,18 

128,126 

128,128 

table-1 


,  this  are  others  programs  with  Matlab  (with  your  graphic  for  analyze  the  behavior). 

Next  to  this,  generates  the  matrix,  we  generates  other  matrix  in  3-d  and  put  the  anterior  matrix  , 
how  the  floor  to  this  one. 

Now  the  procedure  is  calculate  with  the  numerical  method,  iteration  method,  the  values  to  the 
voltages  in  the  other  points  in  the  net,  different  to  zeros  . 

Here  we  have  positions  for  to  solve  in  the  borders  and  in  the  interior  to  the  ‘imaginary’  cube. 


The  method  work  graphically  in  this  form: 


d2V  d2V  d2V 

-2  +  +  ^2  -  o 

^  X  d  y  o  z 

Ax  =  Ay  =  Az  =  Ah 


(W_ 

& 


a 


Vj  -Vo 
Ah  ^ 


F4-Fo 

Ah 


&  , 


V5-Vo 

Ah 


In  similar  manner  make  it  the  other  derivates,  point’s(d,f,c,b  and  e), 
next  have: 


d2V 

I  &  ( 

'dV\ 

dX2 

«?xJ 

V'-Vo-V,^/ 

7m 2 


v ;  +v3-  2V0  / 

/(Ah)2 


with  the  similar  procedure,  we  calculate  the  other  two  axis  (y,z),  and  have  in  the  final: 
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Vt  +  V3-2  VO  +  V4  +V2 -2V0  +V5  +V6-2Vq  =  0 


Vt  +V3  +  V4  +V2  +V5  +V6  -  6Vq  =  0 


Vi 


+v2  +v3  +v4  +v5  +v6/ 

/6 


We  make  the  programs  relative  a  this  problem  and  I  find  a  “solution”,  the  iterations  or  the 
number  of  times,  has  many  variations,  the  variation  range  are  between  [10-100]  loops,  we 
obtain  various  graphics  showing  the  behavior  to  the  simulation.  We  know  that  the  behavior  is 
good,  analyzing  the  graphics.  The  betters  behaviors  was  analyzed  with  100  iterations,  yours  can 
see  the  next  graphic. 


We  look  in  this  graphics,  that  have  it  a  smooth  slope,  the  peaks  correspond  to  values  generated 
by  the  “randomly  peak”  function  ,  the  range  to  the  function  are  between  [+10,-10]  volts,  and  the 
extremes  values  in  the  anterior  picture  are  +10  and  -9.24  volts. 

Other  question  in  this  part  are  to  maintain  fixed  the  points  showed  in  table-1;  Only  the  points 
around  every  these  anterior  points,  will  be  modified  .  With  this  anterior  process,  we  will  fill  the 
floor  to  the  cube,  every  128x128  points  will  has  a  value.  The  graphic  to  these  fixes  points  are  : 
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Now  we  work  with  the  mother  matrix  (mvol)  and  will  fill  it  in  the  “space"  with  16  levels,  graphing 
the  first,  eight  and  last  levels,  the  first  level  is  the  same  that  the  anterior  “smooth  peak  matrix” 
(mvol3d).  With  a  program,  I  calculate  the  voltage  in  every  point,  in  this  case,  could  change  the 
values  in  the  “specials  points”,  and  in  the  formula,  use  information  about  the  upper  space,  we 
will  neglected  by  be  zeros  in  all  the  process;  in  this  program  the  principal  formula  is: 

v(i,j,k)  =  (y(i,j,k  -  1)  +  v(/  - 1  ,j,k)  +  v(/  + 1  j\k)  +  v(i,j  -  \,k)  +  v{ij  +  1,/r))  /  5  . 

this  formula  has  changes  if  we  work  in  the  corners  or  in  the  borders  to  the  spaces,  modifying  the 
numbers  of  factors  by  3  or  4  (dependent).  In  this  program  we  recalculates  many  times  using 
iterations  in  the  loops  hoping  by  betters  solutions  and  more  exact  values  in  every  point. 

The  results,  was  good,  we  maintain  the  smooth  slope  and  the  figures  are  corrected  in  every  layer 
with  the  iterations  without  error  in  the  borders  (we  has  many  in  the  process),  our  last  proof  was 
with  40  iterations  and  finding  the  information  in  the  16  layers,  the  three  graphics  show  the 
advantage  using  these  mathematical  method,  and  the  figures,  repeat  again,  look  very  good. 
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Next  to  the  anterior  process,  we  add  noise  to  the  signal  and  repeat  a  similar  process  that  the 
anterior.  We  work  with  a  noise  signal  generated  with  this  equation: 

noise  =  (2.5*((2*rand(16))-1)), 

when  rand()  is  the  Matlab  random  function,  [0,1],  the  random  data’s  interval. 

Now  the  final  signal  =  “pure"  signal  +  noise  signal 

The  “pure”  signal  is  generated  with  a  smooth  peak  function,  in  this  simulation. 

The  process  was  similar  but  the  result  are  different,  logically. 

in  the  reality,  this  is  the  kind  of  signal  with  that  we  work  in  the  future,  remember,  all  this  process 
is  only  computer’s  simulation’s,  when  will  go  to  the  experimental  table,  or  to  the  space,  we  will 
work  with  many  types  to  noise  and  with  other  types  of  turbulence. 

With  this  process  have  this  pictures: 


and  the  other  graphic,  with  a  comparison  images  picture, 
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Now  we  make  the  process  with  the  16  layers  and  the  graphics  to  the  second,  eight  and  sixteen 
layers  are  in  the  following  pictures  . 

Yours  can  analyze  the  images  process  between  the  first  layer ,  the  eight  layer  and  the  sixteen 
layer, 


control  graphic  z=8  Iterations^  40 


analyzing  the  last  10  rows 


control  graphic  z=16  lterations=  40 


0  0 


yours  can  see  the  change  and  how  convert  these  method  a  distortion  graphic  in  a  smooth  and 
clear  picture,  the  mathematical  method  was  functioning. 

Now  the  process  will  consist  in  to  construct  the  “electronic  part”,  we  will  go  now  to  the  Real 
world,  I  need  now  to  construct  the  hardware  for  all  the  process,  to  make  all  these  simulations, 
“real”. 

I  will  construct  two  parallels  boards,  one  square  128x128  and  other  hex127,  these  boards  are 
immersed  into  a  solution  of  cupric  sulfate(CuS04). 

In  the  principle,  the  boards  are  fixed  in  the  same  axis,  but  next  can  it  moved  one  into  the  space, 
the  system  has  movement  in  the  X  axis  with  one  plaque  fixed  and  the  other  moved  (hex127 
plaque),  moving  in  the  Y  axis  (up,  down),  one  fixed  again  and  the  other  with  movement  and  for 
last  in  the  Z  axis  (left-right),  how  a  car  drive,  again,  one  fixed  and  the  other  with  movement. 

I  need  too  a  mechanical  system  for  control  the  movement  to  the  board. 

Yours  can  see  the  next  schematic  picture,  and  in  there  explain  to  yours  graphically,  how  work  in 
the  real  world  the  summer  project,  in  this  case  put  the  in’s  and  out’s  in  the  intermediate 
positions  to  the  schema  . 
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The  picture  is  this: 


For  this  process  ,  the  solution  between  the  boards  (CuS04,  cupric  sulfate),  will  server  to  us  how 
medium  to  the  “  electric  field  movement"  ,  if  when  we  put  near  or  fair  the  boards,  this  solution 
will  work  how  medium  and  iterpolator. 

I  will  construct  the  two  terminals  (one  for  transmit  other  receiving),  one  for  every  board,  the 
“electricity”  go  out  by  the  128x128  terminal  and  will  receive  it  the  hex127  terminal,  and  the 
electric  waves  “surfer”  into  the  cupric  sulfate,  when  the  process  are  finish,  will  see  the  changes(if 
have  it)  to  the  electric  field  and  we  can  find  the  “perfect  distance”  between  the  boards  and  we 
know  ,  what  values  has  between  the  different  separations,  and  will  reconstruct  the  initial  signal 
(square  128x128)  in  a  final  signal  (hex127),  we  need  too  construct  others  boards,  one  for  square 
and  other  for  the  hexagesimal  board,  with  this  boards,  we  make  the  amplification  to  the  signal 
and  transform  a  initial  signal  to  16  in  a  256  signal,  we  will  transmit  the  “pulses"  and  next  to  the 
capture,  we  make  other  board  and  make  the  remap  in  a  hexagesimal  127. 

How  chemical  information,  we  have  that  the  volume  susceptibility  to  the  water  is 
-90x1  O'6  and  for  the  cupric  sulfate  is  CUS04*5(H20)  =+176x1  O'6 

This  information  will  need  it  to  calculate  the  perfect  distance  and  the  electric  field  make  it  how  a 
constant  value. 

Others  chemical  properties  are: 

1 .  The  cupric  sulfate  is  soluble. 
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2.  The  cupric  sulfate  is  soluble  into  the  water 

3.  The  chemical  reaction  by  which  it  is  formed  is  called  a  precipitation  reaction. 

4.  The  electrical  properties  to  the  elements  are: 

The  last  are  only  a  few  properties  to  the  solution. 

Re-writing  about  the  electronic  part,  we  can  say  that  will  need  op-amp  in  the  boards,  in  the 
moments  don’t  have  a  electronic  schema  to  the  boards. 
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Conclusions 


Maybe  the  conclusions  about  my  summer  job  are  only  a  few  words,  maybe  the  best  or  the 

betters  conclusions  are  when  I  finish  all  the  project,  including  the  electronic  part  and  makes  the 

finals  and  “real"  measures,  the  following  punctual  information  reflex  my  personal  conclusions  to 

this  part  to  the  project: 

1.  Sometimes  the  mathematical  method,  in  this  case,  iteration  method,  could  help  to  the 
solution  for  our  physics  problems. 

2.  In  this  case,  we  can  look  the  process  graphically  and  we  could  saw  how  the  figure  change  in 
every  layer,  to  a  not  very  smooth  initial  graphic,  transform  in  a  very  smooth  and  “clean”  final 
graph(with  and  without  noise). 

3.  If  you  compare  the  noises  graphics  and  the  without  noise  graphic  you  can  say  that  the  initial 
graphics  are  very  ,  very  different  but  in  the  final  had  “similar”  smooth  pictures,  saw  the  right 
bar  in  the  graphic  and  saw  to  you  the  ‘VOLTAGES”  intervals  ,  are  both  are  similar. 

4.  In  the  real  simulation,  we  will  hope  that  the  cupric  sulfate  work  how  the  medium  into  the 
boards,  and  the  “layers”  in  the  computational  simulations,  are  the  different  distances  in  the 
hardware  world. 

5.  We  can  reconstruct  initials  signal,  wavefront,  by  example,  in  pulses,  and  reconstruct  in  any 

form,  in  this  case  a  16x16  128x128  in  16  different  layers. 

6.  Matlab  are  a  very  good  software-packet  for  solve  mathematical/physics  problems. 

7.  The  team  work  is  very  important  to  solve  any  problem,  I  received  here,  many  and  useful 
help,  thank  you  for  all. 
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Abstract 

Magnetized  target  fusion  provides  a  method  to  reach  fusion  ignition  conditions  that  is  intermediate 
to  inertially  confined  and  magnetically  confined  fusion  schemes.  Although  a  magnetic  field  is  used 
to  inhibit  thermal  conduction  from  the  plasma  to  the  confining  liner,  plasma- wall  mixing  still  poses 
a  problem.  When  this  mixing  occurs,  increased  impurities  in  the  plasma  reduce  the  plasma  tem¬ 
perature  and  prevent  the  plasma  from  reaching  ignition  conditions.  Computer  simulations  studying 
the  effect  of  plasma-wall  mixing  due  to  the  Rayleigh-Taylor  instability  in  an  imploding  cylindrical 
linear  are  presented.  The  simulations  are  performed  using  the  2D  magnetohydrodynamic  code, 
MACH2. 
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I.  Introduction 

Achieving  controlled  nuclear  fusion  has  been  a  goal  of  the  scientific  community  for  several  decades.  Commonly 
employed  schemes  include  Inertial  Confinement  Fusion  (ICF)  and  magnetic  confinement  (MFE).  Typical  ICF  methods 
compress  a  gas  to  sufficient  density  and  temperature  to  initiate  fusion.  For  the  minimum  laser  energy,  the  implosion 
process  must  occur  on  time-scales  shorter  than  the  time  required  for  thermal  conduction  into  the  capsule  wall.  This 
requires  expensive  driving  systems  which  deliver  the  necessary  energy  on  nano-second  time-scales.  Compressing  a 
non-magnetized  plasma  on  micro-second  time-scales  provides  time  for  the  high-temperature  plasma  to  transfer  energy 
to  the  environment.  This  reduces  the  plasma  temperature  and  hence,  the  fusion  rate.  Yet  the  presence  of  a  sufficient 
magnetic  field  within  the  plasma  can  reduce  the  thermal  conduction  between  the  plasma  and  its  surroundings.  This  is 
the  idea  behind  MFE  designs.  An  alternative  approach,  whose  density  and  time  scales  fall  between  those  of  ICF  and 
mfe,  is  Magnetized  Target  Fusion  (mtf). 

Forming  a  union  between  ICF  and  MFE,  MTF  is  a  two-step  process.  The  first  stage  requires  the  formation  of  a  warm 
(-100  eV),  magnetized  (-100  kG)  wall-confined  "target”  plasma  prior  to  implosion.  The  "target"  plasma  is  then 
quasi-adiabatically  compressed  by  an  imploding  liner  during  the  second  stage.  The  joining  of  these  two  methods 
requires  mating  a  plasma  formation  system  with  a  target  implosion  driver.  At  the  Air  Force  Research  Laboratory— 
Phillips  Research  Site,  a  plasma  formation  system  has  recently  been  combined  with  a  quasi-spherical  liner  implosion 
system  [1].  Since  the  magnetized  plasma  reduces  the  rate  of  thermal  losses,  electrical  pulsed-power  technology 
(which  operates  on  micro-second  time-scales)  can  be  used.  Such  pulsed-power  implosion  experiments,  similar  to 
those  required  for  MTF,  have  already  been  conducted  [2-4].  Thus,  the  necessary  technology  and  experimental  know¬ 
how  for  MTF  exists. 

The  principle  feasibility  issue  for  the  MTF  concept  is  the  mixing  of  the  fuel  and  boundary  material.  When  a  plasma  is 
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in  contact  with  a  material  surface,  heavy  ions  can  be  sputtered  off  the  material.  These  sputtered  ions  are  highly  ionized 
and  are  responsible  for  large  amounts  of  radiation  losses.  This  increased  radiative  loss  cools  the  plasma  and  reduces 
the  chance  of  reaching  ignition  conditions.  Although  the  magnetic  field  is  intended  to  provide  a  barrier  to  plasma 
particles  reaching  the  liner,  instabilities  can  develop  between  the  plasma  and  the  liner  that  reduce  the  effectiveness  of 
the  magnetic  field.  A  well-known  source  of  non-uniform  behavior  in  electromagnetic  implosions  is  the  Rayleigh- 
Taylor  instability.  In  ordinary  hydrodynamics,  this  instability  occurs  when  a  heavy  fluid  is  supported,  against  the 
force  of  gravity,  atop  a  light  fluid.  The  interface  between  the  two  fluids  becomes  rippled  and  allows  the  heavy  fluid 
to  mix  into  the  light  fluid.  In  MHD,  the  roles  of  the  heavy  and  light  fluids  can  be  assumed  by  the  plasma  and  magnetic 
field,  respectively.  The  Rayleigh-Taylor  instability  occurs  when  the  acceleration  points  from  the  heavy  fluid  (plasma) 
to  the  light  fluid  (magnetic  field).  During  implosion  processes,  the  interface  between  the  plasma  and  the  liner  can 
become  unstable  between  the  times  of  maximum  inward  liner  velocity  and  maximum  compression.  In  the  simulations 
considered  here,  there  is  no  magnetic  field  in  the  plasma.  For  this  case,  the  plasma  plays  the  role  of  the  light  fluid 
and  the  aluminum  liner  that  of  the  heavy  fluid.  The  present  work  considers  plasma-wall  mixing  in  a  cylindrical  liner 
due  to  the  Rayleigh-Taylor  instability. 


Using  the  2-1/2  dimensional  magnetohydrodynamic  code,  mach2  (Multiblock  Arbitrary  Coordinate  Hydrodynamics, 
2D)  developed  at  the  Air  Force  Research  Laboratory  -  Phillips  Research  Site,  a  cylindrical  liner  imploding  on  a 
deuterium  plasma  is  simulated.  Section  II  provides  a  brief  description  of  mach2.  Section  III  details  the  simulation 
model.  The  simulation  results  are  presented  in  Section  IV.  Conclusions  and  future  simulation  considerations  are 
given  in  Section  V.  The  input  deck  used  for  one  of  the  simulations  is  provided  in  the  appendix. 
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II.  MACH2  Description 

mach2  [5]  is  a  finite-difference  code  that  solves  the  time-dependent,  single-fluid,  multi-temperature  MHD  equations. 
Using  an  Arbitrary-Lagrangian-Eulerian  multi-grid  scheme,  the  problem  geometry  is  constructed  from  multiple 
blocks  (logically  rectangular  collections  of  cells)  of  arbitrarily  shaped  quadrilateral  cells,  in  either  planar  or  axi- 
symmetric  coordinates.  The  code  calculates  all  three  components  of  the  velocity  and  magnetic  field,  though  no  vari¬ 
ation  is  allowed  normal  to  the  problem  geometry  (i.e.,  no  z-dependence  in  planar  problems  or  theta-dependence  in 
cylindrical  problems).  The  finite-differencing  is  carried  out  according  to  a  finite-volume  technique,  along  with 
second-order  Van  Leer  convection  and  flux-conserving,  constrained  transport  for  magnetic  advection  [6].  Solving  the 
resistive  MHD  continuity,  momentum,  energy  and  magnetic  field  equations,  mach2  closes  the  system  of  equations 
through  the  equation  of  state,  of  which  several  options  are  available.  Among  these  are  analytic  and  tabular  models 
(sesame  tables  at  Los  Alamos  National  Laboratory),  which  include  the  transport  coefficients.  Several  options  are 
available  for  radiation  modelling.  Multiple  circuit  solvers  are  also  available. 


III.  Simulation  Model 

The  simulation  model  [7]  considers  a  30  cm  high  aluminum  liner  with  an  inner  radius  of  4.9  cm  and  wall  thickness  of 
0.1  cm.  The  liner  encloses  a  deuterium  plasma  with  a  particle  density  of  1017  cm'3'.  The  liner  and  plasma  originally 

r 

have  a  temperature  of  27  meV  and  300  eV,  respectively.  The  liner  is  connected  to  an  RLC  model  of  the  SHIVA  star 
capacitor  bank.  At  t=0,  the  applied  voltage  is  84  kV  with  an  inductance  of  39  nH  and  resistance  of  1  m£2.  The 
deuterium  properties  were  determined  from  the  SESAME  tables,  except  for  the  resistivity  and  thermal  conductivity  for 
which  a  Spitzer  model  was  used.  For  the  aluminum,  tabular  values  were  used  for  the  resistivity  and  thermal  conduc¬ 
tion  coefficients,  along  with  a  Steinberg-Guinan  elastic-plastic  model.  For  simulations  run  in  Eulerian  mode,  a 
hydrogen  gas  was  allowed  to  flow  through  the  outer  boundary  (i.e.,  radius)  of  the  grid.  This  prevents  the  volume  left 
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behind  by  the  moving  liner  to  be  filled  with  a  non-physical  aluminum  gas.  The  hydrogen  gas  employed  a  tabular 
equation  of  state  with  constant  electric  diffusivity  (104  m2/s)  and  thermal  conduction  coefficients  (104  W/m-eV).  The 
gas  was  brought  in  at  0.2  eV  and  0.01  kg/m3. 


Preliminary  simulations  were  run  in  a  Lagrangian  fashion  to 
allow  the  grid  to  follow  the  liner.  These  simulations  used  the 
initial  conditions,  as  described  above,  and  ran  for  25  pis.  The 
problem  geometry  was  scaled  from  30  cm  down  to  1  cm;  the 
plasma  and  the  aluminum  were  each  dimensioned  with  8  cells 
in  the  radial  direction  and  16  in  the  axial  direction.  As  indi¬ 
cated  in  Figure  1,  the  liner  reached  a  2:1  compression  ratio  and 
a  maximum  inward  velocity  at  ~18  pis.  At  approximately  22.5 
pis,  maximum  compression  occurred  with  a  radius  of  about  1.5 


time  (ps) 


Figure  1  Radial  momentum  (in  a  1 -radian  slice  of  the  full  271, 
axi-symmetric  simulation)  from  a  preliminary  study.  The  temperature, 
velocity  and  mass  density  profiles,  as  well  as  the  magnetic  field  and 
current,  at  t  =  18  ps  were  used  as  initial  conditions  for  the  stability 
simulations.  For  the  run  shown,  the  maximum  inward  velocity 
occurred  at  t  =  18.1  ps  and  the  maximum  compression  at  t  =  23  ps. 


To  reduce  computation  time,  the  stability  simulations  were  started  at  18  (is,  using  the  results  of  the  preliminary  runs 
as  "initial  conditions".  The  temperature,  mass  density  and  velocity  profiles,  as  well  as  the  magnetic  field  and  current, 
at  18  (is  were  used  to  initialize  the  problem.  These  were  run  for  an  additional  5  us  (the  case  of  64  axial-cells  was  run 
for  4  ps)  in  an  Eulerian  fashion  to  prevent  "twisting"  of  the  grid.  In  order  to  stimulate  instability,  a  sinusoidal  mass 
density  perturbation  was  added  throughout  the  volume  of  the  aluminum.  In  each  case,  the  perturbation  amplitude  was 
±10%  with  four  complete  cycles  in  the  axial  direction. 

Clear  indicators  of  mixing  between  the  plasma  and  aluminum  melted  off  the  liner  are  the  radiation  cooling  rate  and  the 
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neutron  production  rate.  An  emission  radiation  model  was  used.  It  computes  the  time  rate  of  change  of  the  specific 
internal  energy  according  to  the  opacity.  The  aluminum  was  allowed  to  radiate  only  for  cells  in  which  the  aluminum 
density  was  less  than  100  kg/cm3.  The  neutron  production  model  considers  only  those  cells  that  are  predominately 
deuterium.  For  each  such  cell,  the  plasma  mass  density  and  temperature  are  used  to  compute  the  number  of  neutrons 
produced.  This  model  is  not  considered  to  be  extremely  realistic,  but  does  provide  qualitatively  acceptable  results. 

For  the  Rayleigh-Taylor  instability,  the  fastest  growing  modes  are  those  with  the  shortest  wavelength.  Thus,  the 
cell-size  governs  which  modes  will  be  present  in  the  simulation.  Since  the  proposed  MTF  geometry  has  a  height  of  30 
cm,  studying  a  mode  with  a  125pm  wavelength  would  require  at  least  240  cells  in  the  axial  direction,  as  well  as 
sufficient  cells  in  the  radial  direction.  Problems  with  such  a  large  number  of  cells  require  extremely  long  simulation 
run-times.  In  order  to  the  reduce  the  number  of  cells  and  shorten  the  execution  time,  the  height  was  scaled  to  much 
shorter  lengths.  This  change  in  height  reduces  the  liner  inductance  and  ultimately  alters  the  current  flowing  through 
the  liner.  In  order  to  maintain  the  same  current  as  would  flow  through  a  30  cm  liner,  an  option  was  added  to  MACH2 
that  scales  the  inductance  and  the  voltage  along  the  circuit  path.  This  option  is  called  htsclfac  and  occurs  in  the 
CODtrl  namelist  (It  is  applied  only  to  the  first  circuit  element.).  Defaulted  to  1,  htsclfac  need  only  be  set  if  scaling 
is  desired.  Using  a  htsclfac=T0  means  that  a  30  cm  problem  can  be  simulated  using  a  simulation  height  of  3  cm. 
This  reduction  in  the  total  number  of  required  cells  for  a  particular  resolution  allowed  significant  improvements  in 
execution  times. 


IV.  Simulation  Results 

To  discover  when  the  instability  becomes  significant,  three  different  axial  cell  dimensions  were  considered:  32,  64 
and  128  cells  per  centimeter.  Figure  2  shows  the  radial  momentum  for  each  run,  including  a  run  with  no  perturbation 
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and  32  cells  per  centimeter,  and  indicates  that  the  pertur¬ 
bation  delayed  the  time  at  which  the  liner  pinched  for  each 
case.  However,  a  significant  instability  was  not  obtained 
until  a  mode  wavenumber  of  64ft  cm"1  was  considered 
(i.e.,  128  cells  per  centimeter).  The  reason  for  the  lack  of 
compression  can  be  discerned  from  the  radiation  cooling 
rate,  shown  in  Figure  3.  At  ~20.5ps,  spikes  develop  in  the 
radiation  caused  by  aluminum  ions  sputtering  off  the  liner 
and  mixing  with  the  plasma.  These  ions  then  radiate  sig¬ 
nificant  amounts  of  energy  and  cool  the  deuterium.  Corn- 


Figure  2  Radial  momentum  for  various  Rayleigh-Taylor  instability 
modes.  For  k  =  32n  and  below,  the  perturbation  has  only  a  minor  effect  as 
indicated  by  the  case  when  no  perturbation  was  applied.  For  the  k  =  64n 
case,  the  instability  is  highly  evident.  The  "bump”  at  ~21.5p$  indicates  a 
drop  in  plasma  pressure  due  to  plasma-wall  mixing  lowering  the  plasma 
temperature.  The  delayed  "pinch-time"  resulted  from  a  reduced  liner 
inductance  which  decreased  the  driving  current.  The  larger  momentum 
values  at  1 8  ps  for  the  perturbed  cases  is  because  the  perturbation  added 
more  mass  to  the  liner  compared  to  the  non-perturbed  case. 


Figure  3  Base- 10  logarithm  of  the  scaled  radiation  cooling  rate.  The 
radiation  rate  is  scaled  by  multiplying  the  rate  by  the  ratio  of  problem 
height  to  0.5  cm  (the  k  =  16n  and  non-perturbed  height);  this  gives  each 
case  the  same  volume  properties.  The  perturbation  has  essentially  no 
effect  on  the  radiation  cooling  rate  until  the  k  =  64n  case.  Then, 
considerable  amounts  of  energy  are  lost.  The  reason  for  the  larger 
radiation  rate  for  the  k  =  32 n  at  the  start  and  end  of  the  simulation  is 
uncertain.  (Wavenumbers  are  given  in  terms  of  1/cm.) 


Figure  4  Base-10  logarithm  of  the  scaled  neutron  production  rate.  Tht 
scaled  rate  means  that  the  production  rate  is  multiplied  by  the  ratio  of 
problem  height  to  0.5  cm  (the  k  =  16ir  and  non-perturbed  height);  this 
gives  each  case  the  same  volume  properties.  At  -21.5  ps,  there  is  a 
sudden  drop  in  neutron  production  for  the  64 n  case  due  to  increased 
radiative  losses.  The  larger  production  rates  for  I6n  and  32rt,  compared  to 
the  unperturbed  rate,  are  due  to  larger  temperatures  prompted  by  the 
perturbation;  this  is  probably  a  result  of  the  method  by  which  the  neutron 
production  rate  is  calculated.  (Wavenumbers  are  given  in  terms  of  1/cm.) 


paring  the  radiation  cooling  rate  and  the  neutron  production  rate  (see  Fig.  4)  shows  that  the  plasma  has  lost  sufficient 


energy  that  the  neutron  production  drops  several  orders  of  magnitude.  This  drop  in  plasma  temperature  results  in  a 


loss  of  plasma  pressure.  With  a  lower  pressure,  the  liner  experiences  less  resistance  to  compress  the  plasma  which 
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causes  the  change  in  concavity  of  the  radial  momentum  at  21.ps,  shown  in  Figure  2.  The  liner  then  pinches  much 
later  than  22.ps.  This  is  related  to  the  loss  of  material  from  the  liner,  reducing  the  liner  inductance  and  lowering  the 
driving  current.  With  a  lower  current,  the  liner  is  driven  more  slowly  and  thus  takes  longer  to  pinch. 


V.  Conclusion 

Magnetized  target  fusion  (mtf)  is  a  promising  possibility  for  break-even  fusion  production.  Although  the  presented 
results  lack  the  primary  appeal  of  MTF  designs  (i.e.,  magnetized  plasma),  they  are  necessary  for  an  adequate  under¬ 
standing  for  the  effects  of  high-Z  contaminants  in  the  plasma;  such  plasma  pollutants  can  occur  even  in  a  magnetized 
target  plasma.  The  simulations  demonstrate  that  a  major  source  of  instability  during  the  compression  phase  is  short- 
wavelength  perturbations.  These  can  result  from  nonuniform  mass  density  or  uneven  current  flow  within  the  alumi¬ 
num  liner  .  Furthermore,  the  complexity  of  the  simulations  and  the  numerous  physical  processes  considered  show  that 
mach2  can  successfully  simulate  complicated  phenomena  with  confidence. 

The  next  step  in  the  MTF  modelling  process  is  to  magnetize  the  plasma.  MACH2  already  has  a  built-in  option  (FRC- 
frese)  to  supply  the  desired  magnetic  configuration.  Another  consideration  would  be  employing  mach2’s  material 
tracking  interface.  This  would  allow  the  liner  to  be  followed  in  detail  during  the  development  of  an  instability.  Also, 
development  of  a  more  discerning  neutron  production  rate  routine  would  allow  more  quantitative  result,*  to  be 
obtained. 
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Appendix 

This  appendix  contains  the  MACH2  input  deck  for  the  k  =  64ti  simulation  presented  above.  It  was  run  using  a  modified 
version  9805.  The  modifications  include  a  pwlinear option  for  roinit,  uinit,  vinit  and  teinit. 


Scaled  MTF.  Multigrid.  Strength.  Eulerian.  Perturbation. 
pert02 
Jcontrl 
t  =  I8.0e-6, 
dt  =  l.e-15, 
dtmax  =  1.0e-10, 
twfn  =  23.0e-6, 
cyl*  1.. 
eoson  =  .true., 

(split  =  0, 
radiate  =  .true., 
radmodl  =  ’emission’, 
ciron  =  .true., 
neutrons  =  .true., 
thmldif  =  .true., 
anisot  =  .true., 
mgmodet  =  ’converge’, 
nthrmax  =  2000, 
flxlmt  =  0.4, 
tdtol  =  l.e-4, 
tdrelax  =  1.5, 
hdiff  =  .true., 
mg  mode  =  ’converge’, 
itmaxrd  =  1000, 
ndtol  =  1  .e-4, 
rdrelax  =  1.5, 
aresfdg  =  0.05, 
aresvac  =  l.e4, 
hallon  =  .false., 
magon  =  .true., 
brbzon  =  .false., 
it  pot  =  10, 
potrelx  =  0.25, 
hydron  =  .true., 
theb  =  1., 
volratm  =  0.8, 
courmux  =  5.0, 
rmvolnn  =  0.9, 
itopt  =  20, 
mu  =  5.6, 
eps  =  l.e-6, 


strength  =  .true., 
meshon  =  .true., 
n  smooth  =  4, 
wrelax  =  0.25, 
nigen  =  0, 
niter  =  3, 
eqvol  =1., 
zeroghcl  =  .false., 
multgrd  =  .false., 
mglmax  =  4, 
trackon  =  .true., 
htsclfac  =  240.0, 

Send 
Soutput 
dtrst  =  2.25e-6, 

!  display  every  500 
intty  =  ’edits,50’, 
intty(2)  =  *0’, 
ncyctty  =  500, 
dtp  =  0.1e-6, 
plot(0)  =  ’cpt’, 
plot(l)  =  ’grid’, 
pIot(2)  =  ’velocity’, 
plottype(2)  =  ’vector’, 
plot(3)  =  ’velocity*, 
plottype(3)  =  ’contour’, 
plotcomp(3)  =  ’vectorz’, 
plot(4)  =  ’ni’, 
plottype(4)  =  ’contour’, 
plotcomp(4)  =  ’scalar’, 
plot(10)  =  ’radheat’, 
plonype(lO)  =  ’contour’, 
plotcomp(10)  =  ’scalar’, 
plot(15)=  ’density’, 
plottype(15)  =  'contour’, 
plotcomp(15)  =  ’scalar’, 
p!ot(16)=  ’pressure’, 
plony  pe(  16)=’  contour’ , 
plotcomp(16)  =  ’scalar’, 
plot(17)=  ’thrmflux’ 
plottype(17)  =  ’vector’. 
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plot(18)  =  ’material’, 
plottype(  18)  =  ’contour', 
plotcomp{  1 8)  =  ’scalar’, 
plot(19)=  ’tcone’, 
plottype(19)  =  ’contour', 
plotcomp{19)  =  ’scalar’, 
plot(20)  =  ’nrate’, 
plottype(  19)  =  'contour', 
plotcomp(19)  =  ’  scalar’, 
plot(21)  =  'sigfac', 
plottype(19)  =  ’contour’, 
plotcomp(  19)  =  ’scalar’, 
kpltreg(l)  =  1,2, 
kcon(l)=  11, 
intbound  =  .false., 
clabpict  -  .true., 
fichfram  =  1, 
dtslic  =  0.1e-6, 
sltce(O)  =  *cpt’, 
slice(15)  =  'tcone', 
slice(16)  =  ’nrate’, 
sliced  7)  =  ‘sigfac’, 

Iblkslic  =  1, 
ibdyslic  =  4, 
ijslic  =  8, 
ncychist=  100, 
histenrg  =  .true., 

Send 

Scumt 

circtype(l)  =  ’rlc’, 

!  values  at  t  =  1 8  us 
current(l)  =  4.654e6, 
volt(l)  =  -2.877e4, 
exind(l)  =  8.497e-8, 
exres(l)  *  l.e-3. 

Send 

Sezgeom 

!  height  was  30cm 
nblk  =  2, 
npnts  -  8, 

pointx(l)  =  0.0e-2,  pointy(  1)  =  0.000e-2, 
pointx(2)  =  2.1e-2,  pointy (2)  -  O.OOOe-2, 
pointx(3)  =  2.5e-2,  pointy(3)  =  O.OOOe-2, 
pointx(4)  =  0.0e-2,  pointy(4)  =  0. 125e-2, 
pointx(5)  =  2.  le-2,  pointy(5)  =  0. 125e-2, 
poimx(6)  =  2.5e-2,  pointy(6)  =  0. 125e-2, 
comers(  1 , 1 )  =  4,5,2, 1 , 
comers(I,2)  =  5,6,3,2, 

Send 
Sezphys 
matnamig  =  ’d’, 
siecoldg  =  -l.e99, 
resmodlg  =  ’tabular’, 
etamaxg  =  l.e4, 

Send 
Smatmdl 
matnam(l)  =  ’d’, 

!  ...except  use  hydrogen’s  category  3  table:  resistivity,  etc. 
mattabs(3,2)  =  25252, 


!  ...and  in  this  case  use  spitzer  resistivity. 
resmodl(l)  =  ’spitzer’, 
tcnmodl(l)  =  ’spitzer’, 

!  setup  A1 

matnam(2)  =  ’al-new’, 
resinodl(2)  =  ’tabular’, 
tcnmodl(2)  =  ’tabular’, 
elpmodl(2)  =  ’steinb-g’, 
elpname(2)  -  'al-1 100*, 
rocrad(2)  =  1.0e2, 

!  setup  inflow  gas 
matnam(3)  =  ’gas’, 
sesanam(3)  =  'h', 
eosmodl(3)  =  ’tabular’, 
tcnmodl(3)  =  ’constant’, 
etc0(3)  =s  l.e4, 
resmodl(3)  =  ’constant’, 
eta0(3)  =  l.e4. 

Send 

Sbfield 

momfld  =  2.  le-2. 

Send 
Simnesh 
icells(l)  =  42, 
jcells(l)  =s  16, 
roi(l)  =  7.02e-3, 
tempi(l)  *  965.2, 
teinit(l)  =  ’pwlinear’, 
teir(l,l)=  0.0,0.01725,0.019125, 
teiprfKU)=  1.0,0.9998,  0.998, 
ui(l)  =  -3.415e3, 
uinit(l)  =  ’pwlinear’, 
uir(l,l)  =  0.0,0.021, 
uiprfKl.O^O.O,  1.0, 
magxybc(l.l)  =  ’conmutv’, 
magxybc(3,l)  =  'contnutv’, 
gdvl(l)  =  0.0, 
tethrm(l,l)  =  0.1, 
rothrm(l,l)  =  8.8e3, 
gdvlb(2,l)  =  0.0, 
tethrai(3,l)  =  0.1, 
rothrm(3,l)  =  8.8e3, 
hydbc(4,l)=  ’axis’, 


icells(2)=  8, 
jcells(2)=  16, 
matnami(2)  -  ’al-new’, 
roi(2)  =  2.7e3, 
roinit(2)  =  ’pwlinear’, 

roird.2)  =  2.  le-2, 2.15e-2, 2.2e-2, 2.3e-2, 2.4e-2, 2.5e-2, 
roiprfr(I,2)  =  0.28,  0.65,  1.00,  0.99,  0.97,  0.95, 
tempi(2)  -  4.0e-2, 
teinit(2)  =  ’pwlinear’, 

teiK  1 ,2)  =  0.02 1 ,0.022,0.0225,0.023,0.0235,0.024,0.0245,0.025. 

teiprfi<l,2)=  1.26,1.00,1.07,  1.24,1.47,  1.72,2.00,  2.15, 
ui(2)  =  -3.523e3, 
uinit(2)  =  ’pwlinear’, 

uir(  1 ,2)  =  2.0625e-2, 2.125e-2, 2.4375e-2, 2.5e-2, 
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uiprfr(l,2)=  1.0,  0.979,  0.920,  0.914, 

magxybc(l,2)  =  ’contnutv’, 
magxybc(3,2)  =  ’contnutv’, 
gdvl(2)  =  0.0, 
gdvlb(2,2)  =  0.0, 
hydbc(2,2)  =  ’flowthru’ 
matbc(2,2)  =  ’gas’, 
roflow(2,2)  =  l.e-2, 
tflow(2,2)  =  0.2, 
magzbc(2,2)  =  ’insulatr’, 
currcir(2,2)  -  1, 
gdvlb(4,2)  =  0.0, 
ropert(2)  =  ’roderick’, 
roprtam(2)  =0.10, 
roprtmy(2)  =  4.0, 
roprtln(2)  =  0.125e-2, 
binit(2)  =  ’annulijy’, 
bzi(2)  =  109.0, 

Send 


I 
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Abstract 

Spatial  light  modulators  are  used  in  a  wide  variety  of  applications  ranging  from  optical  computing  to 
display.  The  application  that  is  of  interest  for  this  report  is  that  of  real-time  holography.  A  system  architecture  has 
been  developed  that  uses  a  spatial  light  modulator  as  a  real-time  holographic  element  in  a  telescope  system  to 
perform  adaptive  optical  correction  of  primary  mirror  distortions.  This  allows  a  less  expensive,  lightweight  primary 
mirror  to  be  used  with  little  or  no  degradation  of  the  resulting  image.  The  spatial  light  modulator  is  the  key  element 
to  this  system.  A  phase  conjugate  projection  system  is  also  being  investigated  for  beam  collimation/projection 
applications  using  a  similar  inexpensive,  lightweight  primary  mirror. 
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IdgSfiPPS  aberration  compensation 

Modem  telescope  designs  require  the  use  of  large  primary  mirrors  to  achieve  the  high-resolution  images 
desired.  In  order  to  maintain  the  high  optical  quality  necessary  for  these  mirrors,  they  are  generally  manufactured  as 
a  heavy,  monolithic  component.  Fabrication  of  such  a  component  is  very  costly,  and  transporting  it  is  very  difficult. 
Because  of  these  factors,  the  size  of  ground-  and  space-based  telescopes  is  limited.  By  using  an  adaptive  optics 
technique,  the  optical  requirements  for  the  primary  mirror  can  be  relaxed  so  a  large,  low  optical  quality  primary 
mirror  can  be  used  in  conjunction  with  some  compensation  optics  to  achieve  high-resolution  images  with  a 
considerably  reduced  cost. 

One  proposed  method  of  adaptive  optic  compensation  for  low  quality  primary  mirrors  is  by  using  real-time 
holography  with  an  optically  addressed  spatial  light  modulator.  This  technique  uses  holographic  phase  subtraction 
to  remove  the  aberration  of  the  low  quality  mirror.  Figure  1  shows  a  conceptual  set-up  of  the  holographic  phase 
subtraction  technique.  A  beacon  laser  probes  the  mirror  and  picks  up  the  aberration  of  the  primary  mirror,  <t>a,  where 


with  X  is  the  wavelength  of  the  beacon  laser  and  OPD  is  the  optical  path  difference.  The 


aberrated  beacon,  Be^a ,  is  interfered  with  a  reference  plane  wave,  B ,  from  the  beacon  laser  to  form  a  hologram. 
For  an  amplitude  hologram  the  transmission  is  proportional  to  the  squared  modulus  of  the  illuminating  field,  given 

by: 


T  ~  2  +  ei4>a  +  e~'^ 


This  hologram  is  then  illuminated  by  the  aberrated  object  field,  Ae’^e"^  where  Ae,<l5  is  the  undistorted 
object  field.  Multiplying  the  aberrated  object  beam  with  the  transmission  function  of  the  hologram: 

Ae*V*aT~2ei<V,|>a  +  Ae'^e2'^  +Aeio 
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lightweight  mirror  with 
phase  aberration,  e1*3 


Figure  1.  Schematic  of  an  adaptive  optical  technique  for  phase  subtraction. 

The  three  terms  represent  the  zero,  plus  one,  and  minus  one  diffracted  orders,  respectively.  The  zero  order, 
proportional  to  el0e^a ,  is  the  original  aberrated  object  beam  transmitted  through  the  hologram.  The  plus  one  order, 
proportional  to  eia>e2l<^a ,  is  the  addition  of  the  aberrated  phases,  creating  further  distortion.  The  minus  one  order, 
proportional  to  el0 ,  is  the  subtraction  of  the  aberrated  phases,  leaving  only  the  unaberrated  phase  information,  O  . 


Real-time  holographic  element 

The  key  component  to  this  whole  procedure  is  the  real-time  holographic  element,  the  OASLM.  There  are 
several  important  characteristics  to  consider  when  thinking  about  dynamic  compensation  of  larger  aberrations  and 
severe  distortion,  namely  sensitivity,  speed,  diffraction  efficiency,  and  resolution.  OASLM  devices  have  been  used 
to  compensate  for  large  distortions  (hundreds  of  waves  of  aberration)  with  beam  intensities  as  low  as  50  pW/cm2. 

Aberration  compensation  requires  that  the  OASLM  react  at  a  speed  greater  than  that  at  which  the  aberration 
is  changing,  which  is  very  dependent  on  the  specific  application.  Nematic  liquid  crystal  limits  the  response  time  to 
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about  30  Hz  due  to  the  slow  recovery  time  but  devices  have  been  built  that  have  high  diffraction  efficiency,  -30%. 
Surface  stabilized  ferroelectric  LC  can  have  a  much  faster  response  time,  >lkHz,  but  typically  have  much  lower 
diffraction  efficiency,  <10%,  because  of  the  low  tilt  angle  of  the  molecules. 

Since  the  compensated  image  is  found  in  the  first  order  diffraction,  the  diffraction  efficiency  of  the 
OASLM  is  also  a  major  source  of  loss.  For  a  phase  only  sinusoidal  modulation  the  refractive  index  profile  can  be 
written  as  n  =  n  +  An  cos  qx .  The  diffraction  efficiency  of  this  index  profile  is: 

r)  =  |J](27tAnd/A.)|2 

where  27iAnd /X  is  the  phase  associated  with  the  transmittance,  or  half  the  total  phase  modulation.  With  a  nematic 
OASLM  operating  in  the  analog  phase  modulation  mode  the  largest  theoretical  diffraction  efficiency  achievable  is 
only  about  34%  at  a  phase  modulation  of  about  3.7.  In  addition,  phase  modulation  occurs  for  only  one  polarization 
direction  in  nematic  liquid  crystal,  so  for  unpolarized  readout  light  the  diffraction  efficiency  is  reduced  to  about 
1 5%.  If  a  binary  phase  grating  could  be  written  in  a  polarization  independent  material,  the  diffraction  efficiency 
could  theoretically  go  as  high  as  40% 

The  resolution  of  the  OASLM  is  limited  by  the  diffusion  of  charges  in  the  photoconductor  and  by  electric 
field  fringing  effects  in  the  photosensor,  reflecting  layer,  and  liquid  crystal  layer.  The  resolution  of  the  OASLM 
limits  the  severity  of  distortion  that  can  be  corrected  for.  It  is  desirable  to  have  a  resolution  of  greater  than  100  line 
pairs  per  mm  while  maintaining  high  diffraction  efficiency,  however,  commercially  available  OASLMs  have  greatly 
reduced  diffraction  efficiencies  at  spatial  resolutions  above  about  30  lp/mm. 

A  prototype  device,  which  has  been  built  under  contract  to  the  US  Air  Force  by  the  Research  Institute  for 
Laser  Physics  and  Peterlab  in  St.  Petersburg,  Russia,  uses  a  new  type  of  photoconductor  and  a  deformed  helix 
ferroelectric  liquid  crystal  (DHFLC)  in  order  to  improve  the  resolution,  diffraction  efficiency  and  the  response  time. 
A  schematic  of  this  device  is  shown  in  figure  2.  The  photoconductor  is  carbon-doped  amorphous  hydrogenated 
silicon  (a-Si:C:H)  which  has  a  very  low  dark  conductivity  as  well  as  low  carrier  diffusion  .  The  liquid  crystal  layer 
is  a  ferroelectric  liquid  crystal  (smectic  C*)  with  a  high  tilt  angle  of  40°  aligned  in  deformed  helix  configuration. 

The  photoconductor  allows  an  increase  in  the  resolution  while  the  DHFLC  increases  the  speed  and  diffraction 
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Figure  2.  Schematic  depiction  of  the  prototype  optically  addressed  spatial  light 
modulator. 


efficiency.  In  addition,  if  the  high  tilt  angle  of  the  ferroelectric  liquid  crystal  makes  it  less  sensitive  to  the 
polarization  of  the  read-out  beam. 


This  device  can  operate  as  a  binary  phase  modulator  with  the  application  of  a  driving  voltage  and  the 
writing  light.  Initially  the  DHFLC  is  driven  to  one  of  its  stable  states,  ±40°,  by  the  application  of  a  large  voltage 
(initialize  voltage)  across  the  device,  which  results  in  a  large  electric  field  across  the  liquid  crystal  layer.  The 
polarity  of  this  voltage  is  arbitrary,  the  effect  is  the  same  regardless  of  the  starting  state.  Then  the  voltage  is 
switched  to  a  small,  opposite  value  and  the  writing  light  is  allowed  to  strike  the  photoconductor.  Where  light  hits 
the  photoconductor  the  local  conductivity  increases  and  the  electric  field  across  it  drops,  therefore  the  electric  field 
across  the  liquid  crystal  layer  increases  and  causes  the  molecules  to  switch  to  the  other  stable  state  rapidly.  This 
process  is  illustrated  in  figure  3.  The  areas  not  illuminated  by  the  writing  light  also  switch  to  the  other  stable  state, 
albeit  on  a  slower  time  scale,  causing  a  decay  of  the  written  pattern.  This  decay  can  be  avoided  if  the  entire  process 
is  repeated  before  the  molecules  in  the  dark  regions  have  time  to  completely  switch  states.  The  response  time  of  the 
device  is  the  time  it  takes  to  reach  maximum  diffraction  efficiency.  The  refresh  rate  of  the  OASLM  is  determined 
by  the  repetition  rate  of  the  device,  not  by  the  response  time  of  the  liquid  crystal  itself. 
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Figure  3.  Timing  diagram  for  hologram  formation  in  the  prototype  OASLM 
showing  the  bias  voltage,  spatial  and  temporal  profile  of  the  writing  intensity 
pulse  and  the  orientation  of  the  index  ellipsoid. 

Ideally  the  light  regions  would  switch  rapidly  and  the  dark  regions  would  maintain  their  orientation  for  an 
extended  time  creating  a  persistent  hologram.  Physically,  however,  there  is  a  tradeoff  between  the  response  time 
and  the  hologram  persistence.  In  order  to  have  a  fast  response  time  there  must  be  a  large  electric  field  across  the 
liquid  crystal.  So  faster  response  times  require  higher  applied  voltages  that  also  result  in  higher  electric  fields  in  the 
dark  regions  and  a  faster  decay  rate. 


Results 

A  primary  focus  of  the  investigation  is  the  speed  and  diffraction  efficiency  of  the  OASLM  when  writing 
and  reading  a  hologram.  Figure  4  shows  the  experimental  set-up  used  to  characterize  these  properties.  A  doubled 
Nd:YAG  (532  nm)  laser  was  used  as  the  writing  beam  and  a  diode  laser  (810  nm)  was  used  as  the  readout  beam. 
The  two  different  wavelengths  were  used  because  this  device  doesn’t  have  a  reflective  layer  and  must  be  used  in  the 
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Figure  4.  Experimental  configuration  used  for  characterizing  the  prototype 
OASLM 


transmission  geometry.  In  order  to  avoid  the  readout  beam  being  absorbed  in  the  photoconductor  and  degrading 
performance,  a  higher  wavelength  was  used.  The  writing  beam  was  passed  through  a  liquid  crystal  shutter  to  control 
its  temporal  profile  and  then  defocused  and  sent  through  a  Mach-Zehnder  type  interferometer.  The  two  arms  of  the 
interferometer  were  overlapped  on  the  photoconductor,  creating  a  straight-line  interference  pattern  and  phase  grating 
in  the  liquid  crystal.  The  readout  beam  was  diffracted  from  this  grating  and  the  temporal  dynamics  of  the  first  order 
were  measured  using  a  fast  photodiode.  Because  the  device  has  no  anti-reflection  coating  and  there  was  still  some 
absorption  of  the  read  beam  in  the  photoconductor,  the  diffraction  efficiency  is  defined  as  the  power  in  the  first 
diffracted  order  divided  by  the  power  in  the  transmitted  zero  order  when  there  is  no  grating  present. 

The  resolution  of  the  device  is  measured  by  changing  the  crossing  angle  of  the  writing  beam  and  therefore 
the  spatial  frequency  of  the  grating  being  written  and  observing  the  resulting  diffraction.  The  diffraction  efficiency 
can  also  be  changed  by  varying  the  values  of  the  applied  voltages  and  the  parameters  of  the  writing  pulse.  Figure  5 
shows  a  family  of  curves  that  were  obtained  by  changing  the  resolution  of  the  pattern  being  written  on  the  OASLM. 
For  these  curves,  the  resolution  was  set  at  350  line  pairs  /  mm  and  the  voltage  and  write  pulse  were  adjusted  to 
achieve  the  maximum  peak  diffraction  efficiency.  These  settings  were  then  left  constant  and  the  resolution  was 
changed  by  adjusting  the  crossing  angle  of  the  writing  beams.  The  solid  lines  represent  the  diffraction  efficiency, 
the  write  pulse  is  the  dotted  line  and  the  voltage  is  the  dashed  line.  The  decay  in  efficiency  was  discussed  earlier 
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Figure  5.  Transient  diffraction  efficiency  as  the  grating  resolution  is 
changed. 

(c.f.  fig.  3)  and  the  sharp  drop-off  of  efficiency  is  due  to  the  voltage  switching  back  to  the  large  initializing  value.  It 
has  been  shown  that  by  changing  the  amplitude,  duty-cycle  and  frequency  of  the  driving  voltage  it  is  possible  to 
make  the  diffraction  efficiency  nearly  constant  over  the  entire  read-out  time. 

Figure  6  shows  the  maximum  peak  diffraction  efficiency  as  a  function  of  the  spatial  frequency  of  the 
grating.  The  maximum  peak  efficiency  was  found  by  making  small  adjustments  to  the  voltage  and  write  pulse  at 
each  spatial  frequency.  A  peak  diffraction  efficiency  of  34%  can  be  attained  at  18  lp/mm  and  at  370  lp/mm  a 
diffraction  efficiency  of  8%  is  still  achieved.  Note  that  the  resolutions  were  calculated  from  the  angle  between  the 
zero  order  and  first  order  diffracted  beams  ■ 
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Figure  6.  Maximum  peak  diffraction  efficiency  as  the  spatial  grating 
constant  is  changed. 


Phase  Conjugate  Projection 

It  is  sometimes  desirable  to  have  a  large  scale  beam  collimator  for  use  with  a  telescope  system  for  such 
uses  as  illuminating  an  object.  Figure  7  shows  a  schematic  of  a  system  that  uses  phase  conjugation  with  a  large,  low 
optical  quality  projection  mirror  to  achieve  beam  collimation.  As  in  the  system  above,  the  aberrations  of  the 

primary  mirror  are  probed  with  a  beacon  laser  resulting  in  a  wavefront  of  Be1^3 .  The  aberrated  probe  beam  is  then 

directed  into  a  degenerate  four- wave  mixing  phase  conjugate  mirror.  The  phase  conjugate  of  the  probe  beam,  with  a 

•/ 

wavefront  of  Be-1^  ,  is  projected  onto  the  primary  mirror  and  the  aberrations  in  the  conjugate  beam  cancel  the 
aberrations  on  the  mirror  leaving  only  a  plane  wave.  At  this  time  the  investigation  of  this  system  is  just  beginning 
and  there  are  no  publishable  results. 
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Conclusion 

For  applications  involving  real-time  compensation  of  large  dynamic  aberrations  requires  the  use  of  a  spatial 
light  modulator  beyond  the  capabilities  of  what  is  currently  available  commercially.  A  prototype  ferroelectic  SLM 
was  built  under  contract  to  Phillips  Laboratory  to  meet  these  capabilities.  This  deformed  helix  ferroelectric  liquid 
crystal  SLM  shows  increased  diffraction  efficiency  and  a  much  larger  resolution  than  existing  devices.  For  lower 
resolutions  the  diffraction  efficiency  can  be  made  as  high  as  33%.  In  addition  to  having  high  diffraction  efficiency, 
the  prototype  device  shows  a  very  high  resolution,  at  350  line  pairs  per  mm  the  diffraction  efficiency  was  still  near 
10%!  This  device  shows  characteristics  that  would  make  it  a  very  good  candidate  to  be  used  to  compensate  for 
severe  aberrations. 
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